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ABSTRACT 


This thesis contains discussion, theory and program code 
mora computational fluid dynamics (CFD) model of a wing of 
Pebitrary planform. The model assumes incompressible, 
/mviscid, ireotational flow. The program computes forces 


meting on the wing by modeling the” flow with a set of horse 


shoe vortex elements. It models the flow over an arbitrary 
forme using two solutYrons. One solvcton “isthe ideal lift, 
associated with a cambered and twisted wing. The other 


solution is the additional lift associated with a flat wing. 
The program computes wing camber and twist using an elliptic 
loading distribution. The thesis includes the FORTRAN source 
code, a separable User’s Manual for the VORTEX program, 
esccussionm of the theory applied in the model, and 
imvetructions for operating the program. It shows a sample 
wing planform with tabular and graphic results. 

The thesis also discusses two other CFD models based on 
circulation (I) and pressure difference (ACD). It presents 


some of the problems and solutions in grid generation. 
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I. INTRODUCTION, PURPOSE AND GOALS 





S20 > Deersaranm introductory course in aerodynamics for 


aeronautical engineering students at Naval Postgraduate 


Sehool, @enterey, CA. The course must provide students with 
the required knowledge fre futwre courses ieee t he 
Aeronautical Engineering Curriculum. The student must 


understand the techniques used in computational fluid 
auoamres (CPD) to satisfy that requirement. Powerful 
computational models are available jee g the working 
mem@odynamicist, but they are seldom usable iy ea 
imme roductory “course. The programming language for these 
models is usually optimized for computational efficiency 
rather than teaching, and the documentation may be poor or 
unavailable. Students need a computer program with graphic 
output which is less complicated to use as a teaching tool. 
Providing that computer program is the goal of the thesis. 
ive plroreram funs On a micro or desk top computer, in 
keeping with the emphasis on use of individual workstations 
to supplement the school’s mainframe computer. The “arte © o 
computer has better graphic output programs available, 
without the extensive programming required on the school 
Hemmitame. This flexibility allows the student to modify the 


output to suit his/her requirements. 


The thesis contains five principal sections. The firs 
section, TECHNIQUES, provides a brief synopsis of three 
different flow models. The section also provides’ some 
insight to certain problems that arise ineGroe 

The second section, VORTEX PROGRAM DEVELOPMENT, contains 
a detailed discussion of the horse shoe vortex model. The 
section includes the theory and techniques used in the 
VORTEX program. 

The jhe section, CIRCULATION MODEL, contains 
discussion of a circulation (1) model The section includes 
the theory behind the model and some of the problems 
encountered. 

The fourth section, PRESSURE DIFFERENCE MODEL, contains 
discussion of the pressure difference (AC) model. The 
section also includes the theory behind the model fama 
problems encountered. 

The fifth and last principal section, the USERS MANUAL, 
is separable from the body of the thesis and provides 
detailed theory on the horse shoe vortex model. The section 
also includes instruction on use of the program andi 
example wing planform with the resulting data. 

The computer listing for the VORTEX program=eis tiewauaer 


ln the-Appendix,. 


A. 


PROGRAM ASSUMPTIONS AND REQUIREMENTS 


To keep the model streamlined, the following assumptions 


were made and requirements established. 


ite 


Flow is steady, inv LS curd , irProtational, ati Ga 
incompressible. 

The wing surface has zero thickness. 

The planform has any sweep, taper, and aspect ratio. 
The program must support low aspect ratio planforms 
Gress ™tnan i270). 

The program must contain comments, for easy 
modification and change. 

The program results must be available in tabular and 
graphic form. 

The program must use accepted variable names. (aspect 
melo, “caper ratio, etc. ) 

The program must have two options. One is to generate 
loading for a flat wing shape. Mite: Notier sas: ato 
generate a wing shape for an elliptic loading. 

The programming language must be FORTRAN, with 


executable and source files provided. 


Il. TECHNIQUES 


Three different models were investigated. Each of Witeme 
models uses a governing equation for the induced velocity at 
each. “conturod point! on the wing. That velocity (ime 


associated with the vortex strengths located at field points 


on the wing. The induced velocity is a function of Mga 
factors. One is the vortex strength and the other Toueene 
physical orientation of the field and controle pounce Flow 


must be tangent to the wing surface and is found by adding 
the induced velocity vector (w) to the remote velocity 
VECCOr (Ve). The end result is a direct correlation between 
the wing's shape? and the strengths of the vortices. This 
relationship allows either the strengths of the vortices or 


the wing shape to be the independent variables. 


IField points are the positions of the vortex elementer 
and control points are the positions where the velocity 
evaluated. For example, a control point is where flow 
tangency is enforced. 


2Wing shape will be used throughout to describe the 
suxftace ‘slope of the “wine. For example, one wing shape 


might be a flat plate, another might have the tip twisted 
relative to the root lcno nae 


The three models use emreulatwon (Tae incremental 
elreulation (Al), and pressure difference coefficient (AC)). 


These quantities are related as follows. 








e& (=| 
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ao HORSE SHOE VORTEX MODEL 

This program models the flow by a set of horse shoe 
MmomelLces distributed over the wing. Each vortex consists of 
a bound portion perpendicular to the remote velocity, plus 
evor tCrailing portions. Each vortex element generates an 
incremental circulation (AI) around the element. 

The program generates a matrix of influence coefficients 
from an equation for the induced velocity. The program 
finds the strength of each bound vortex and therefore the 
myeremental circulation around the element by solving the 
set of simultaneous equations for wing shape. The program 
can also invert the problem and solve for the wing shape 
meem a specified vortex distribution. 

The program can then find the forces and moments on the 
poe, Knowing the distribution of incremental circulation. 
haes solution is possible by using the Kutta-Joukowski 
theorem which relates incremental circulation and effective 


velocity to the force generated. 


B. CIRCULATION MODEL 

This program models the flow by a series of circulation 
elements distributed over the wing. There is also a 
circulation sheet that extends from the trailing edge to 
i mye eye 

Solution of a set of governing equations provides the 
strength of each circulation element. The equations satisfy 
flow tangency on the wing surface and two boundary 
COnNGLeloms. Firstly, circulation strength is Zero “alone 
the leading edge and wing tips. Secondly, the derivative of 
circulation with respect to chord-wise direction is zero at 
the trailing edge. The derivative of circulation wire 
FESPeCe (66 CROnd mice olen) ty. This trailing edge boundary 
condition satisfies the Kutta condition of Zero vores 
strength at the trailing medce- Zero vorticity at” Veue 
trailing edge also means that the circulation strength 


behind the wing is a function of the span-wise coordinate 


Only: 


C. PRESSURE DISTRIBUTION MODEL 

The final model of the flow uses a series of pressure 
difference elements, ACy, distributed over the wing. The 
program builds a set of N equations in N unknowns from a 


governing equation and from the requirement for flow 


tangency. The boundary conditions are: 
i ACy is zero at the trailing edge, and along the wing 
tips. 


2. For a flat wing, the strength of AC) aon 2 thre 
leading edge is unrestrained. 
Pesoen cumin = of wwe llipure Tift distribution, the 


strength of ACy along the leading edge is zero. 


D:. PROBLEMS INHERENT TO ALL MODELS 
sie Field and Control Point Placement 


The placement of field and control points within 


elements is important for all models. Goimeident Ereld and 
control points exhibit singularities. Each model uses 
different methods to handle them. Direct. Integration 
eliminates the singularity in two models. The other 


requires a special arrangement of control and field points 
along the leading edge and wing tips. 


Dini Grid Size and Solution Speed 


Fine grids require more time for solution and larger 


computers. All models compromise between data quantity and 
Speed. The circulation model is very memory intensive and 
will only run on the mainframe computer. The other models 


will run on a desktop computer, though the time required 
increases rapidly with finer grids. 
3. Grid Element Shapes 

Finite elements in a non-rectangular planform can be 
either trapezoids, or rectangles. Trapezoids match the 
planform exactly; rectangles present a ragged leading and/or 
Grailing edge. The form to use is dependent on the type of 
model. The circulation model can use a trapezoid without 


7 


di iti cmee,.. Coordinate transformations easily handle the 


shape and at the same time allow non-uniform cosine spacing 


of the elements. The Vortex and ACy models work best with 
rectangular grid elements. Therefore the planform is not an 
exact representation of the wing. This modification jae 


reasonable since the forces and moments computed are 


comparable to other models. 


E. SCOPE OF DISCUSSION 

Of the three models investigated, only the VORTEX 
program was completely successful. Significant effort wae 
expended on the CIRCULATION and AC, models trying to make 
them functional. Rather than expend a large portion of this 
thesis discussing all aspects of the two unsuccessful 
models, an overview of them will be presented and the 
relevant factors of each will be addressed. Since the 
VORTEX program was successful, and is completely functional, 
a detailed discussion of the procedure and assumptions is 
provided. In keeping with this, the FORTRAN source code for 


the VORTEX model is the only computer code included. 


III. VORTEX PROGRAM DEVELOPMENT 


A. BASIS OF THE VORTEX MODEL 

The VORTEX program is an adaptation of a program by 
Moran [{Ref.1]. The program models the flow by a series of 
horse shoe vortices distributed over the wing. Icne ti i s 
text, Moran [Ref.1] develops a simple program to find the 
strengths of a series of horse shoe vortices. His program 
is rather limited in that it only works for straight wings 
without taper or sweep. raddteronatiy his program only 
Finds the loading generated by a flat wing. The final 
Iemitation is his use of uniform sized grid elements, which 
fails to concentrate grid elements near the boundaries of 
the wing. 

As with many aerodynamic problems, the VORTEX program 
uses a grid or mesh of finite sized elements distributed 
over the surface of the wing. The program divides the wing 
planform (there is zero thickness modeled in this program) 
into N discrete elements. 

The program solves a set of N equations for N unknowns. 
The equations are derived from an induced velocity equation 
and the requirement for flow tangency on each grid element. 
The N unknowns are the strengths of the individual horse 


emoe vortices associated with each element. 


Knowing the strengths of the vortex elemen cicero 
program can find the forces and moments acting 6n Chews 


aS a result o£ Ghosewe ements. 


B. MODIFICATIONS TO THE MODEL 

The VORLAT program, developed by Moran [Ref.1], required 
changes to suit the goals of the thesis. Items requiring 
improvement were the ability to handle sweep and taper. 
Other improvements were also added to generate the grid 
Using Cosine (spacing. A module to determine camber and 


twist, called shape, was also incorporated. 


C. TECHNIQUE FOR A FLAT WING 

This is a short description of the VORTEX “prograue 
Consult Moran [Ref.1] for a description of the ~ VOR 
program, the foundation of the VORTEX progran. 

1. Down-wash and the Relationship to Flow Tangency 

A horse Shee vortex is located in the xy plane as 

shown in Figure 3-1. The vortex induces a velocity at any 
point in the xy plane. That velocity can be determined from 


equation 3.1) sasyiollows:. 


eee aa 
ar V (x-x PP +y-y¥)" Ar | V(x-x te ye 


w (2) = anty—y,) 


———{1+ 
aN ye (x—x ) (ca) 
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Rie oaiemx. then 


ay l a 3°29 
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wij (x,y) HomGimem ve TOCLtY indweed ae the control point (x,y). 
ArT is the incremental circulation!?, BUG (ka eae aud «{ Xs 7 f ) 
are the coordinates of the corners of the horse shoe vortex. 


Motan derived this equation using the Biot-Savart Law 


[Ref.1]}. — 





A Horse Shoe Vortex in the Wing Plane 
Bigcure 3-1 





lnifferent authors use different variable names for 


Circulation and Vorticity. In the Vortex model, AY is used 
mere incremental circulation. faethe Girculation model, F is 
Eeeca for Circulation strength. Caution is urged when 


relating one model to another. 


IL 


As shown in Figure 3-2, the remote velocity, with aevame 


at angle of attack a, has components V,cos(a) and Vosi7 (age 


Vo 
Veo sin@®) 





Veg COSRH FLAT WING 


Velocity Components on a Flat Wing 
Fipure 3-2 


If the flow is tangent to the wing, the down-wash 


2 


Vebociey at a point on the wing must equal Vgsin(a). Each 


vortex on the wing contributes to the down-wash”? at every 
point on the wing. So, an equation can be developed for each 
control point as a function of Che Gotrengen sce aach horse 
shoe vortex. In this thesis, the mid point of any bound 
Vortex is termed =a Eteld (poumer Any point on the wing 


\ 


where flow tangency is evaluated is termed a control point. 


2Downwash is considered positive when its directions 
downward. 


3The contribution will be positive when the induced 


velocity is downward, or negative when the induced velocity 
is upward. 


eZ 


Be. 


Pomeexample cf the equation for the control point (Xp ,Yp) Sh, 


Vsina— YT w (x y)=0 


a ae a 


Sombaning all=the control points results in a linear set of 


N equations in N unknowns. 


2. Placement of the Horse Shoe Vortex 


Moran 


placement of a 


with one chord-wise element. 


the vortex is 


[Ref.1] 


lumped vortex 


placed 


others have discussed the 
Ona ew oOo] Gimens tonal airfoil 
For the two-dimensional case, 


Guile tema chondmsno int sand ft low 


tangency is evaluated at the three-quarter chord point. The 


two-dimensional reasoning can be applied to a wing of finite 


span, assuming the wing is 


element. 


CIRCULATION 


c/4 


formed with a single chord-wise 





CHORD 


Crrcwlaetronmisttipucton over a Flat Airfoil 


Pag dem 3-4 


le 


For the flat airfoil, the circulation @iscer toe 
is known to be approximately as shown in Figure 3-3. The 
centroid of this circulation distribution is at the quameee 
chord polnee As a result, it is reasonable to lump the 
total circulation at the quarter <hordwipora For these 
conditions, the flow is tangent at the three quarter chord 
Po Lit - This arrangement correctly computes the momene 
coefficient (C,) for the two dimensional f£lacyate ore 

The reason for other authors’ placement of the 
vortex at the quarter chord of each grid element was 
questioned during development of the VORTEX program. The 
other authors used lifting line theory with a single chords 
wise element for their quarter chord vortex placement. The 
VORTEX program does not use a single chord-wise element. 
The wing is divided into a number of individual elements 
which are modeled with a uniform distribution of vorticity 
over each element. The centroid is at the center of the 
element. oe for the arrangement in Figure 3-4) the 
incremental circulation is concentrated at the center of 


each element. 
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VORTICITY 


CHORD 


PPPpLOncimtationeort Vorticitys with Individual Elements 
Figure 3-4 


WLth@menme “Circulation concentrated at the center of 
each element, the flow tangency point is set at one half an 
element downstream of the concentrated circulation nee falls 
at the border between elements, or at the trailing edge of 
the wing. 

\ 


The coefficients affected by the vortex placement 


are Cy and Xcp. In Ref.2, Xep/e for a straight flat wing, 


aspect ratio 2, was determined to be 0.209. With the lumped 
fometex at the mid chord point, the location of Xep/c 
momverped to a value of approximately 0.234. The error 


Mmiagtecates the vortex is too far aft on the element. 
Relocating the lumped vortex element to the quarter chord 
point, and evaluating flow tangency at the three-quarter 
chord point, moved the computed location of Xep/c forward, 
Pemeestill not to the true value of 0.209. 


Ales) 


In Ref.3, Hough looked at optimum grid Wand 'vorce: 
arrangement for rapid convergence. One of the planforms 
was a straight flat wing with an aspect ratio of 2. Hough 
used a conventional vortex lattice arrangement of uniform 
size grid elements with the lumped vortex at the quarter 
chord” porne. For the aspect ratio 2 wing, he found that 
the computed Cy; /a was larger than actual and that the value 
converged as in Figure 3-5. It was expected that the error 
in C;/a could be reduced by decreasing the planform arcagiag 
some factor. He found that convergence of ©; /a, Vortex iia 
Factor Ore and Xeop was improved by insetting the tip 
vortex by some fraction (d) of a grid element. This inset 
improved convergence dramatically when d = 1/4. Though not 
addressed by Moran [Ref.1l], that is the reason for insetting 
the grid: by de] b7a- 

The VORTEX program uses cosine spacing instead of 
uniform spacing so insetting the grid by d = 1/4 was not an 
available option. Instead, the span-wise grid layout is 
developed with an additional row of tip elements which are 
not used in the computation of incremental circulations 
force. This reduces the planform area and the improvement 


in convergence can be seen in Figures 3-5 and @aae 


4Vortex Drag Factor, K is; m*AR*Cpjy 


16 


paareionally =the convergence of Xeop is improved, as seen in 


Baeswre 3-7. 


12.00 


a 
O 
° 


8.00 


\ ——AR = 2, 4 chord vortices, d = 1 element 
—~——AR = 2, 4. chord vortices, d = OQ 
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Percent error in CL/a 


2.00 
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Vortex Drag Factor (K) 
Figure 3-6 
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-~-—-AR = 2, 4. chord vortices, d = 0 


Oe22 





0.20 
0.00 0.01 0.02 G05 0.04 0.05, 0.06, 0.07 
(Number of Chordwise vortices) 
Xep 
Figure 3-7 
3. Satisfaction of the Kutta Condition 


To satisfy the Kutta condition, the vortex strength 
must be zero at the trailing edge. While not explicieie 
required “by “the VORTE tee proc vam. the resultant vortices 
strengths, for a flat wing, approach zero at the tradi 
edge. Therefore, the results seem to imply that flow 
tangency at the trailing edge of a flat wing is a corollary 
of the Kutta ‘condition. 

4. Planform Development 

It is necessary to develop wing geometry from aspect 
ratio, taper ratio and sweep angle. Aspect ratio is the 
span divided by the average chord. Taper ratio is the )tae 
chord divided by the root chord, and sweep angle is the 
angle between the leading edge and the y axis. The y axis 


is perpendicular to the remote velocity and parallel to the 
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earth. The aircraft design courses taught at NPS use these 
variables to define the planforn. Moran [Ref.1] uses a 
mixed) TrOGtechorawequal to 1.Q0%and varies the span to achieve 
different aspect ratios. The VORTEX program adopts this 


approach, which works well. 


Remote Velocity eA 


am 


Planform and Variables used in VORTEX Program 
Figure 3-8 


The subroutine SET85 determines planform variables 
and stores them in part of the eee WING and SECTN for 
further use. The WING and SECTIN matrices also contain final 
mesults of the program. 

9. Grid Development 

Pic Pp ucrmamendeve lOpc ea Senta “to model the vortex 
system after establishing the outline of the planform. The 
grid elements are rectangular. The number of elements must 
be small to keep the time required for solution reasonable. 
The program concentrates grid elements in areas where the 


rate of change is most rapid, or where the values are most 


1 


iImpertant- Those areas are along the boundaries. An 
excellent method for concentrating the points in the areas 
desired while minimizing the total number of points is 
through cosine spacing. The program uses a cosine function 
to distribute span-wise and chord-wise grid points in a non- 
uniform fashion after a suitable coordinate transformation. 
Figure 3-9 contains a sample of the layout. Notice how the 


rectangular elements model the planform. 


Remote Velocity 





Grid Model with Cosine Spacing 
Figure 3-9 


6. Solving the Set of Equations 





The method used to solve the set of N linear 


equations is arbitrary. Moran [{Ref. 1] uses a Gaussian 
algorithm for solution of an N Dy@eN mate The VORTEX 
program uses the same algorithn. Increasing the “toeam 


number of grid elements beyond 100 would require 


modification of the GAUSS subroutine. 
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ve Finding the Forces on the Elements 


ThemeeEta- Jouwrowski wethneotem states that the force 


per unit span acting on an element is. 


MmGiis| equation, Varrf is the local effective velocity at 
the center of the element. Wee is defined in equation 3.5 
eaaerigure 3-10. p is the density and AI is the incremental 
circulation around the element. Titiscescleecwlation (AL) “Ls 
the same AI contained in equation 3.1, and can be termed the 
incremental circulation that occurs over the element. 

mGeoughout this section, certain approximations and 
dimensional simplifications will be made. The first-is the 
small angle approximation, where the sine is approximately 
equal to the angle in oP radians, and the cosine is 
approximately equal to 1.0. This approximation requires 
that angle of attack for the flat wing be small, generally 
less than 10 degrees. The approximation also requires that 
any wing slope on the elliptically loaded wing be small. 


Mints requirement is met by restricting the desired lift 


eoefficient to values less than about 0.5. In the process 
Oeeecetermining coefficients of lift, drag and moment, 
certain dimensional reference quantities arise. These 


quantities are density (p), remote velocity (V.), planform 


area (S), and average chord (c) . In performing dimensional 


Zi) 


analysis an arbitrarily value of unity can be assignedumam 
these terms. 
a. Finding ACy 
The. distribution ce ACyp over the wing is 
desired, so the force on each element must be converted to a 
dimensionless pressure ditterencewweue omc me nee Vf pe 
defined relative to the wing as; 
— > 


V =) asia 
@® 


Ne 
eff #2 (vsina = w) 338 





Components of Effective Velocity (Viz 
Figure 3-10 


It is also significant to note that AY can be weweten mace 


Aan f 3326 


tags 


Pree the bounasportien or the horse shoe vortex is aligned 
with the y axis. iiemmeesUMemOmmecte force "cross product, 


Mmowaclon 3.4, is; 


2 ; , 
— = o|(w— V,sina jar + V cosa ATR 
NA 


Setting p = 1, Vo in and incorporating he small angle 


approximation; 


— 
AF 2 > oe 
a w- ajar? + are 
Ay 

fms force is related to ACy. ACy ms aupscalawetracher than 


a vector quantity, so the program uses the force component 


normal to the wing to compute ACp. That component is; 


AF = AyAP 3.9 


Making the force non-dimensional gives a pressure 


difference coefficient, AC 


Pp ° 
, AF 
AC = oie ae, 
Pp ] 9 
~pV > AxAy 


73 


Finally, after setting p ="l and) oe ACy for an elememne 


is; 





This equation provides the AC, for each element, which can 
be plotted versus chord, for each span-wise section. 
b. Finding the lift, drag and moment 
Lift, drag and moment are also found frome 
vector force on each element. Equation 3.7 showed that the 
force on an element has components normal and tangent to the 
wing. Once more setting p = 1, and Vy = 1, the vecucm 


force on an element is; 


| y ie 
AF = |(w- sina JATAyi + cosa ATAyk ae 


These forces are oriented relative to the wing coordinaes 


system, where i is parallel to the wing and k is 
perpendieular “tothe ‘wine. Lift and drag are normal and 
tangent, respectively, to the remote velocity (Vo). The 


Wing is at an angle of attack («) to that remote velocity. 


Using this information, the program transforms the forececieae 
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An» 


emictoordinate system oriented to the remote velocity. For an 


iorcividual element the lift andwdrag forces are; 


AL = Elemental Lift Force = (cosa) (AyAI’) (cosa) — (w — sina)(AyAP) (sin a) 


AD = Elemental Drag Force = (cosa) (AyAP) (sina) + (w — sina) (AyATP)(cos a) 


Marne aj small angle approximation, and discarding 


higher order terms, Se a, 


AL = Elemental Lift Force = AyAT — waAyAl = AyAT 


AD = Elemental Drag Force = AyAl'a + (w—a)AyAr = wAyAP 


Mees cCotal lift and total drag acting on the wing is a 
Summation of the elemental lifts and drags. Using wing 


eymme cry ; 


L=Toallif=2 SAL 
semt— span OR hee 


D=ToulDrag=2 SY AD 
semi—span 


The moment about the y axis for an individual 


element, with the accepted sign convention, is; 


AM = Elemental Moment = —(normal wing force) (x position of center of element) 


AM = Elemental Moment = (— AyAT) (x on 


7) 


Total moment about the y axis is a summation of 


the elemental moments. Once again, using wing symmetry. 


M = Total Moment = 2 », AM 3 ate 


semi —span 


8. Finding the Coefficients 


Non-dimensionalizing the total lift and drag gives, 








Cc oa ie 
a - v-s 

a - 

3. 

b= D 
= v°s 

De a 

where S = span x average chord = (b x c), p = 1 and V, =m 


Non-dimensionalizing the total moment about’ the y 


axils gives, 


’ M 
om (a ; 3.18 
we) 6 (ee: 
gts 
where c is the average chord. 
The aerodynamic center X,, is important. Xa. 15 tie 


chord-wise position about which the flat wing has zero 


moment, and is found from. 


x _ — tolalmoment 3. ae 
~ total lift 
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This equation assumes that total moment is taken about the y 
fers, Or x=0. When the program finds the moment on the 
Pebipticalty loaded wing, it uses the X,, from the flat wing 


as the reference axis and finds Cyac. 


D. TECHNIQUE FOR SPECIFIED LOADING 
1. Finding Camber and Twist from Loading 

Elliptic loading in a span-wise direction results in 
minimum induced drag, and. “loading ~s1i1s proportional to 
eureeulation. At sufficiently high aspect ratio, a flat wing 
with elliptic area distribution generates elliptic loading 
mat weighs more than a straight wing. Wings with straight 
leading and trailing edges also cost less to manufacture 
Emam elliptic wings. So, the program generates a span-wise 
eeemptrc Litt’ distribution by slightly twisting the wing. 
The direct relationship between wing shape (camber and 
twist) and circulation distribution makes this generation 


possible. 


ZF 


twist 





Camber and Twist on the Elliptically Loadedqmvanmz 
Figure. 3-10 


For a thin wing, elliptic load distribution ingjege 


chord-wise direction requires an approximately parabolic 


camber shape. The program specifies an elliptic chord-wise 
load distribution and then finds the wings’ associated 
shape. 


2. Specifying the Lift Coefficient 
Specifying the form and amplitude of the ideal lift 
distribution over the wing determines the form and amplitude 
of the wing shape. A load distribution of elliptic form is 
imposed in this case and the corresponding amplitude is 
fixed by the desired ideal wing lift coefficient eG 
on Maximum Desired Lift Coefficient 
The small angle approximation breaks down when large 
Lifte-seoerficien cs are specified: To prevent this 


CeCCUrrence: the maximum desired i it coefficient is 


restricted to 0.5. 
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4. Achieving the Desired Lift Coefficient 


The program scales a reference distribution of ACy, 
which is elliptic in form, to achieve the desired total lift 
meetficient for the wing. 

5. Forces, Moments and Coefficients 

Once the ACy distiawubutimom 1S “scaled for the 
elliptically loaded wing, the program finds forces, moments 
and coefficients in the same manner as for the flat wing. 
The only difference is in the coordinate transformation used 
to convert forces on the wing to lift and drag. The 
individual elements on the elliptically loaded wing are no 
longer at a uniform angle to the remote velocity. This 
variation in angle must be taken into account’ when 
performing the coordinate transformation which gives lift 
and drag forces on individual elements. 

The program finds moment coefficient about’ the 


aerodynamic center, Cya,, from, 


Pimece X.. is found from the flat wing, and ¢ is the average 
elioird . CMa Ss tiie se PlEcChnNe moment of the ‘elliptically 
loaded wing about the y axis, and Cy; is the lift 


Soetticient of the elliptically loaded wing. 


(ie | 


E. VALIDATION OF RESULTS 

Two sources are used to check the validity of the VORTEX 
program results. Lifting line theory provides one source of 
predicted lift and drag values for straight high aspeer 
ratio wings. A continuous loading method developed by the 
National Aerospace Laboratory of the Netherlands (NLR), 
presented in Ref.2, is used as the second source. The 
method used by NLR is a very accurate computational method 
and is considered to be exact for this comparison. 

Kuethe and Chow [Ref.4] discuss lifting line theory for 
the case of a flat, untapered, rectangular wing. Using 
lifting line theory and the VORTEX program, lift and drag 
coefficients for identical wings were computed. A range of 
aspect rations from 0.5 through 20 were selected. Lif @isme 
line theory is known to be reliable at higher aspect ratios, 
but tends to lose accuracy at low aspect ratios. 

As the graphic results in Figures 3-ll and ))3=am 
demonstrate, CL/a and CDi/(a)? values obtained from lifting 
line theory and the VORTEX program converge nicely at high 
aspect watios- At lower aspect ratio, lifting line and Wehbe 
VORTEX program show a difference in calculated value. The 
values of CL/a and CDi/(a)4 obtained by NLR [Ref.2] agree 
closely with the VORTEX program results at aspect ratiouge 
The NLR results validate the VORTEX program results @jmom 
Stralencaswings.. Exact data from swept and tapered wings 


were not investigated. 
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IV. CIRCULATION MODEL 


A. BASIS OF THE CIRCULATION MODEL 

This program models the flow with a sheet of distributed 
circulation and a continuous trailing vortex sheet  biehwag 
the wing. Barna [{Ref.5], and also Milne-Thomson [Ref.6, 


pages 171-177], discuss this conceptually simple model. The 


program uses the wing shape, flow tangency, boundary 
conditions, and the Biot-Savart Law to develop a Sciam 
equations. The program then solves the set of equations for 


the unknown circulation strengths at all grid points on the 


wing. 


B. TECHNIQUE FOR SOLUTION OF THE CIRCULATION MODEL 


dees Induced Velocity 


The velocity induced at a point by the distributed 
circulation sheet on the wing and the trailing vortex Sheer 
can be treated as the sum of the velocity induced by each. 


The velocity induced at the control point (x5,y¥o) Dy oases 


point (x,y), located on the wing ise 











r is the distance between the field and control points, and 
Tl is the strength of the circulation at the field pou 
CXS yore Setting the semi-span length equal to unity, and 


a2 


wae) Eleommchonrdmaccommodates Gititerent aspect and taper 


fae © Ss... 


- - Cma- T> 


Remote Velocity pie = 4 
: 4 
i 


Ewamrormeand Vatrlabucsuused im Cireulation Model 
Figure 4-1] 


Circulation (Ir), and pressure difference coefficient 


* 


(ACp) G@eeerelated as follows: 





Mor equation 4.2, ACy is the pressure difference coefficient 
between the upper and lower surfaces, and [Tf is the 


Sereulation. 


a5 


The velocity induced at control point (x8, )o) >) ee 
segment of the trailing vortex sheet attached to the wing at 


CX 3) ease 








+1 r—(x,,—-x ))/or,,,> 
=| ] | me I — Jas Goa 
trail —1 V aur Sao) dy 


r is the distance between the points, and Ty is the strength 
of the circulation at the field point (x,,y) om Ehewecra tee 
edge. 

When the remote velocity is unity, the induced 
velocity is the same as the slope of the wing at the point 


in question, and the total slope is: 











: 0 wing ° trail 
addletonal yy 
Ww 
(= | wing 4.5 
ox V 
O wing - 
and 
Ww 
(=) = trail 4.6 
ax V 
° trail 2 


The program uses these equations to satisfy flow 


tangency on the wing with a known shape and solves for the 
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ietwown Circulation strengths: fe macmeex  TOtCatlon | the 


equations have the following form. 


cles ea aie: aa tal 67 


2: Boundary Conditions 


The requirement for flow tangency at each element on 
eae, wing iS not enough to find the circulation which 
satisfies the Kutta condition; the program must also enforce 
Sertain boundary conditions. 

The flat wing, which produces the additional lift, 
has the following boundary conditions: 

l. © = 0 along the leading edge and the wing tips. 
2. oé©/éx = 0 along the trailing edge. 

fiiemeecubapercally loadeds wing, which produces the 

deal lift, has the following boundary conditions: 
1. YF = 0 along the leading edge, and wing tips. 
2. @é6F/éx = 0 along both the leading and trailing edges. 

Using a finite difference scheme, it is possible to 
approximate the partial derivatives of [Cf (ér/éx, & Ar/dy) 
from the values of [ at neighboring elements and satisfy 
these boundary conditions. The partial derivative of I with 
Mmespect to x, namely (é@r/dx), approaches infinity at the 
leading edge of the flat wing. This singularity makes 
Mmace differencing near the leading edge difficult. 


Cosine spacing is used to minimize that problem and also to 


SS 


Spacing is used to minimize that problem and also to 
concentrate elements near the edges of the wing. Theses 
coordinate in the chord-wise direction can be expressed in 


an alternate form by ¢, where: 


mantel Pee 
c 


$= cos" ~ 


The additional variables are defined in equation 4.24 and 
Figure 4-3), The y coordinate in the span-wise direction 


becomes @, where: 
0 = ws” '(y) 4.9 


After this “transtormacicn =. 


dl alae 
a = ( <sinig) 4.10 
ap ax \2 

From the transformation, it is possible to show that jae 
has a finite value at the leading edge even when @6r/édx is 


int nice 


36 





evesebe weet beoeeeeboewroen become wiosewete wowed oeceeedoesoedocea=-: 
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eeeeaqdececan 





eseteoeure tf eece te tfoewwrecteoeoteewséeerteee{ mere erstteceeeg eres edoeaean=]: 





Gireurartilon vs Pri 
Figure 4-2 


Using Figure 4-2 as a guide, the program estimates 
the value of <@°F/d¢d@ from the values of [ on neighboring 
elements. The following shows the procedure for finding 
@l/dgdg at the leading edge element as a function of [ on 
neighboring elements. Other points are derived in a similar 


jmashion. 


e Q 3 etal 
r= Ab + Ao" + Asd 


Z 


dl 
7a A, + 2A, + 34,6 


memeche first element, Ad 
For the second element, 2) ea 


For the third element, 2 


oF 


which becomes 


0.5 0.25 0.1257( A,Ag r 


Ee 


135 92.20) 3-3 fo as =7F 


3 
2.5 6.25 15.625) A, se | "oa 


and 


A Ad Sense oe 
Od C ae St, 


} 2 3 


s), 
AZO. =a in Gules \ 
yh US ae 200) a 


I Z 3 


Acad. = tie 
Ad ea econn 


focal 
» p 
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4.132 


4.136 


4.13% 


4.15a 


4.15b 


4 .i135¢ 


when 4.16 is combined with 4.15, it becomes; 


ar l 4.17 
(=) = —{.5r + 66677, -—.19T 
dp/, Ad %) Po %5 
imemuncertainty in this approximation is of order Ag. For 


interior elements, 








aie 1 
(=) = | sr 5 4.18 
dd i Ad Fes eT 
and for the trailing edge element, 
; 1 
(2) = + | rse2n + 95651 - 0431, | peed 
dp/, Ad 9 el 
Pamilar relations hold in the span-wise direction. Jie! 
Matrix form: 
ar fn 0 
COW) = = 
op 
ar Tour yaelt 
BN ) ee 
a0 
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3. Sweep and Taper 


The program must model swept and tapered wings. fixe 
makes a series of coordinate transformations to accomplish 
this. The final result is an equation for the slope due to 


the circulation on theswine 











eat ¥—-Y,\csind (av 4.22 
) sino — } = | do dO 
3 ad 3 2 \ao 


r E 
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where 


= (x ey) 
Our es; 7 By,) 
B=(A - co) 


A = tan(A), tangent ofleading edge sweep 


c =rootchord 
m 


t= 1 — laper ratio 





x, = yho 
c= cl — oy) 


0 = +1 onright wing semi—span 


0 = +1on left wing semi—span 


, - Cna-T > 


menore Velocity See al 


a “ 


4. Field and Control Point Placement 
PreclO anGecol@LOole ponies are located at the center 
of the grid elements. When each grid element contains both 
mere ld and control point, the circulation ehraeeeresults. Loom 


solving the equations oscillates wildly and bears no 
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resemblance to the expected solution. After considetape. 
probing, a modified arrangement was discovered which 
generates a reasonable function. Distribution of a "Cities 
of elements containing only field points along the tips and 
leading edge gives a satisfactory solution. The remainder 
of the grid elements contain field and control points. The 


final (configuration sms. 


[| Field and Control Points 


Remote Velocity 


+ Field Points Only 





Modified Grid Geometry 
Figure 4-3 


This arrangement gives up some degrees of freedom, but that 
can be compensated by adding an extra row of grid elements 
along the chord and an extra column of grid elements along 
the span. 
5. Influence Coefficient Matrices 

Equations 4.22 and 4.23 can be divided into mulewoeee 
integrals. The partial derivatives (dr/d¢g & aP/aa) are 
evaluated at the centers of the elements, and the integrals 
in 4.22 and 4.23 can be represented as a summation over the 
elements. Using these assumptions and equations 4.19 
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fas 


Giows hl 42ers possiblewtomwrite sthesset of equations in 
Matrix form. 

The program builds the final matrix of influence 
Seerricijents from a number of subsidiary matrices. When 
Samation 4,20 and 4.2] are combined with equation 4.22 and 


the cosine coordinate transformation we get. 


Once [M] is generated, it is possible to obtain the 
weaculation solution vector {T). 
6. Solving the Set of Equations 
The program uses the LEQIF subroutine from the IMSL 
FORTRAN library as a linear equation solver. 
7. Finding the Lift, Drag, and Moment 
PeOtm=merc™ Glstribution Of Circulation strength it is 
Pessible to find the lift and drag. Barna {Ref.5] and 
Milne-Thomson [{Ref.6] derive the equations for lift and 
mere, aS a function of circulation and down-wash. The 
Meropram uses these equations to get lift and drag from the 
Seeculation on the trailing edge elements. fie bee Or 2a 


Span-wise section which lies between y and (y+tAy) is, 


AL = pVIAy a a 
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and the induced’ drag on@ay sect Menmarce 


AD = pVI'wAy ee | 


Tr is the circulation along the trailing edge and w is the 
down-wash at the trailing edge. The procedure for f£imdaiae 
moment and center of pressure is not as simple. The program 
would need to resolve forces on the individual elements to 
find moment coefficient and center of pressure. The 
CIRCULATION program was abandoned before Fortran code to 
generate moment coefficient and center of pressure could be 


Wri Cten. 


C. PROBLEMS ENCOUNTERED 
1. Coincident vs Non-coincident Points 

The manner used to handle singularities is one of 
the biggest problems faced when using finite elements to 
model ‘Continuous funcetons:: The circulation model behaves 
very well when the field and control elements are not 
COlMecident. 

Field elements induce velocities which increase Very 
rapidly as the distance to the control point decreases. 
When the field and control points are in the same element it 
appears the induced velocity is indeterminate. However, 
Milne-Thomson [Ref.6] shows that an element induces zero 
velocity (on. (leseire This result of the singulamety oe 
important and requires special @aetencrom 
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Ph Field and Control Points in the Wing Interior 


ihewdsasectionmot velocity induced at a” control point 
Maeitecs) neighboring fdelid poimts is significant. Using a 
simple 3X3 element grid near the center of the wing as an 
example, Figure 4-4 shows the direction of the induced 
welocity from 8 field points. All eight field points have 


merereewvye circulation and surround the center control point. 


REMETE YaeeEOCLLY 
—_— 


4 control point 





[ntrerlor Control Foint 
Figure 4-4 


The three elements upstream from the central control point 


induce a down-wash on the control point. The other five 
elements induce an up-wash on the control point. The net 
velocity is a sum of all eight elements. Recall that the 


velocity induced is dependant on dr/dx, not on T itselfl. 
This net velocity will be down-wash if @ar/dx in the three 
upstream elements is sufficiently greater than in the other 


five elements. 


Ithe derivative of circulation is pPrEOpOELtCIOonal €¢ ACy. 
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3% Field and Gontrol Points ait the Leadingmmarce 


A problem becomes apparent when checking the induced 
velocity at an element on the leading edge. For a 2X3 gaema 
on the leading edge, Figure 4-5 shows the direction of 
velocity induced by each of the neighboring field poimaeuyam 


the central scontrol@poinc 


ReEME TE ae Bee iii 
—_—_—---————————— 


& control point 





Leading Edge Control Point 
Figure 4-5 


For this situation, all the induced velocity is Up -walum 


assuming positive circulation. However, the wing requipecum 
net down-wash for flow tangency. Using this arrangement of 
points, the solution vector oscillates wildly between 
positive and negative values. This oscillation is a Yesumes 


of the system trying to generate down-wash at the leading 
edge elements. 

The program uses a simple solution to this problem. 
By ignoring the requirement for flow tangency on the row of 
elements along the tips and leading edge, the program gives 
a net down-wash at all other elements. This simplificatwem 
results in a well behaved, positive circulation 60Ver seme 


ened retwine: 


It appears the set of equations has more unknowns 
mianm equations, Dweeeooundary, conditions provide the 
meaditional constraints. 

4. Model Complexity 

Developing a computer program usable as a teaching 
tool for basic aerodynamics is the principal goal of the 
thesis. The CIRCULATION program does not fully achieve that 
goal. The concept is relatively understandable, but is 
faeeticult to execute. The complex method needed to build 
the influence coefficient matrix reduces its usefulness as a 
teaching tool. The student needs a simpler tool. The 


eB TEX program discussed in Section III is that tool. 
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V. PRESSURE DIFFERENCE MODEL 


A. BASIS OF THE AC, MODEL 

Like the vortex and circulations modelce the pressure 
difference model (ACy) also uses a governing integral 
equation as its foundation. The integral equation relates 
wing slope at a specific point on the wing to the 
distribution of pressure difference over the whole wing. 
Dividing the wing into N grid elements, with known slopes, 
the program solves for N pressure differences (ACD). 

The requirements of the Kutta condition are iia 
explicitly satisfied by the equation. Instead, additional 


constraints are necessary at the wing tips and trae 


edge. The simplest form of satisfaction is to require that 
ACy be zero at those grid elements. That constraint is more 
Severe than necessary. ACy must be zero at the edge of the 


element but the value at the center of the edge elements can 
be finite and is approximated by fitting a polynomial 
through the edge and neighboring control  podntce This 


approach is similar to that used in the CIRCULATION modem 
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B. TECHNIQUE FOR SOLUTION OF THE ACy MODEL 


5 The Singularity 


The imeciecalemequatton “poverming the wing 


slope/pressure difference relationship is, 


Z (x,y) = | k(x—x,,y—y,) AC, (x,,y,) dx, dy, 5.1 
S 


where the slope, 


° 0 
Z (x,y) = = IZ, (20) 





and, 


] 


Sr -y,)° 


(x—x,) 


k 





|e 








V (x=x,)° + (y—y,)° 


mmen the field and control points are coincident (r=#0) or 
share the Same y value (y=yl), there is a strong 
Paeaeularity. heme composcuple: CO evaluate this integral and 


eliminate the singularity. 
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oe Evaluation of the Integral 


Using the mirror image element on the opposite semi- 


span, it is possible to reduce the integral) tae 


Z Gy) = | | k'(x—x.,y—y,) AC, (x. ,y,) dx, dy 
sn l 1 Plea 1 1 


where, 


k =kt+k 5a 
and 
1 1 (x—x,) 


— Lo Sea6 


8r (yty,) V (z=x,)" + (y+y,) 





The singularity was eliminated by evaluating the integmeame 


Using Figure 5-1 as a guide, the result is. 


2 
REMOTE VELOCITY 
> v1 

e 


+i | 


Variables used for a Grid) 2) veneu. 
Figure 5-1 
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mck”, the mirror image is; 


: ] 
a 3, ot Rio ty, + @,) 


Sa. 


where, 


Qo oe 


qa - Wt ye 
and all other values are the same aston The end resus 
1S: 
kk =k+h° 
and 


Liga = » k” AC, (x,,y,) 
S/2 


In matrix form, the equations can be wrieten: 


He = FI 


A computer program will rapidly calculate all of these 


values. 
Se The Kutta Condition 
Solution of the equations does not guarantee a 
result which satisfies the Kutta condition!. To enforce the 


Kutta condition, the program sets AC, at elements along the 


1 It should be noted that the requirement for flow 
tangency at the elements near the trailing edge was 
sufficient to satisfy the Kutta condition in the )VOHisas 
program. However that result was discovered after the 
PRESSURE program was abandoned and no effort was macdeueee 
extend that reasoning to the PRESSURE program. 


3) 


wing tip and trailing edge (the nth elements) equal to a 
[Ma@etion Ofethe value of the [n-l]th elements. The program 
enforces the condition by fitting a polynomial through the 


two points. See Figure 5-2. 


RoMeis Vv GCHEIT Y 
a 


No eania N- th Elements 
Figure 5-2 


a Grid Spacing 


The program uses a non-uniform grid spacing, based 


Sa a cosine function. The elements near the wing tip are 


small. This spacing has the advantages noted earlier. 


C. PROBLEMS ENCOUNTERED 


In matrix form, the model has the form. 





Peeeeis the influence coefficient matrix, and 2° is the wing 


slope. The model specifies some values of ACy on the left 


ao 


side of the equality and some values of wing slope (2°) on 
the right side of the equality. This complexity prevents 
use of the normal matrix solvers. As a result, the solution 
requires a complex matrix manipulation which is not readily 


understandable, or desirable for a teaching tool. 


D. CONCLUSIONS 

The evaluation of the integral over the element has 
definite advantages for coincident field and control points. 
The evaluation eliminates a strong singularity and should 
give a result using finite elements that is very close to 
the actual property. 

The Kutta condition and method used to £entoiee 


satisfaction is not optimal. Further studies are needed. 
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B. INTRODUCTION 

COMO pio He bM wrote this program, as a tool for use 
ma AE2O35. GieweprOrpram 1S partial satisfaction of the 
requirements for the degree Aeronautical Engineer from the 
Naval Postgraduate School. The objective of the program is 
to provide a simple computer simulation of flow over a thin 
wing with low aspect ratio. High aspect ratio wings are 
better served by a lifting line model. Many advanced 


programs using elaborate methods to model viscous effects, 


compressibility, boundary layer growth, and shocks, are 
available. None of them gives a good introduction to the 
field of computational aerodynamics. This program should 


Fill that need at the Naval Postgraduate School. 

If you wish to skip the theory behind the program, go to 
section I for instruction on running the progran. 

There are a number of ways to model the source of 


Beeces acting on a body surrounded by fluid in motion. 


These include potential functions, vortex distributions, 
Seeceulation distributions, and pressure differential 
mestributions. Each model is related to the other, and each 
has advantages and disadvantages. This program uses a set 


of horse shoe vortices to approximate the flow over a wing. 


5 


C. SYSTEM REQUIREMENTS 


System type IBM PC or compatible 
Memory 256K, minimum 
Math Co-processor (8087) Recommended, Nn Ome 


required. 
Graphtes card Recommended, nl Ome 


wequ hr car 


Graphics program aw Y graphing program 
like GRAPHER, is 
recommended, n Gee 


required. 


D. CONCEPT DESCRIPTION 
l. The Horse Shoe Vortex 
This program, which models the flow by a series of 
horse shoe vortices, distributed over the wing, is an 
adaptation of the VORLAT program by Moran "/ Renee He 
develops a program to find the strengths of a series of 
horse shoe vortices associated with a flat rectangular wing: 
The VORLAT program has its foundation in two-dimensional 
airfoil theory, and Moran [Ref.1] adapted the theory to 
wings of finite aspect ratio. 
The Users Manual will present a short description of 
the VORTEX program. For coverage of the VORLAT program, the 


foundation of the VORTEX program, consult Moran )/pe pa 
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a. Down-wash and the Relationship to Flow Tangency 
The velocity induced at a point in the xy plane 
by a horse shoe vortex, also located in the xy plane, can be 


determined from: 


Be he 2, 2 
Ar V (x—x)) a Oia) | Ar | V (xx) yy) 7.1 


Sa — {)] + RTE rae 
amy —¥) (x—x ) 4n(y—y,) (x— x) 


mo (xX = x,), 


feets the velocity induced at (x,y), called a control point. 
AT is the strength of the horse shoe vortex, and (X3,Y,) and 
(Xa,VYp) are the corners of the horse shoe vortex. The 
center of the horse shoe vortex is a field point. 

The velocity at a point on a flat wing, which 
Meeeat an angle of attack a, has components Vycos(a) and 


Peoin(a) as indicated in Figure 7-1. 


Vea Sind 


Vag COSC FLAT WING 


Velocity Components on a Flat Wing 
Figure /-1l 


ol 


In order to ensure the flow is tangent Cougeme 
wing, the induced velocity or down-wash! at that point Wen 
the wing must also be Vgsin(a). Each vortex on the walmge 
contributes to the down-wash2 at every control point on the 
wing. So, an equation for each control point, as a £unetaee 
of all the field points can Dewvriecen- An example of the 


equation for the control poincuc. yy ise 
F p’Yp 


V_sina — S AP w (x yJ=0 
a ee 


When the control points are combined, the result is a se@tumee 
N equations in N unknowns. 
b. Placement of the Horse Shoe Vortex 
The bound portion of the vortex (that portmom 
perpendicular to the onset flow) is placed at the 1/4 chord 


point of each grid element on the wing. Flow tangency is 


evaluated at the 3/4 chord point of each grid element. 


lnown-wash is considered positive when its direction. 
downward. 


2The contribution will be positive when the induced 


velocity is downward, or negative when the induced velocity 
is upward. 
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c. Satisfaction of the Kutta Condition 
Domestics y secme | Kutta  conmdition, the vortex 
Berenpth must po to zero at the trailing edge. For wings 
and airfoils of zero thickness, flow tangency enforced at 
the trailing edge appears to be a corollary of the Kutta 
Sondition. 
2. Developing the Wing 
a. Planform Geometry 
It is necessary to develop wing geometry from 
aspect ratio, taper ratio and sweep angle. Aspect ratio is 
the span divided by the average chord. Taper ratio is the 
tip chord divided by the root chord, and sweep angle is the 
angle between the leading edge and the y axis. Moran 
{Ref.1] uses a fixed root chord and varies the span as 
necessary to achieve the required aspect ratio and taper 


ratio. The VORTEX program uses this approach, which works 


well. 
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Remote Velocity Bre 


; | 


Planform and Variables used in VORTEX Program 
Figure /-2 


Aspect ratio, AR, is (b2/S). Span (b) isin. 
distance from wing tip to wing tip, and area (S) is the 
total planform area of the wing. 

Taper ratio, A, is the ratio of tip chord to@igoae 


chord. 


tipchord 7a 


— Owes Oe 


rootchord C 


Sweep angle delta, A, is the angle that the leading 


edge makes with a line drawn perpendicular to the root 


chord: 
b. Grid Geometry 
Using aspect ratio, taper “Yateiiog and sweep 
angle, the subroutine SET85 develops a grid. The program 


stores the coordinates of all necessary points in the 
matrices WING and SECTN. Results of the program are also 
stored in WING and SECTN which are written to datageeasiee 


when the program finishes. 
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In a manner analogous to that used by Hough [Ref.3], 
Mae tip wortices are inset to ‘improve accuracy of the 
results. The inset distance is one element wide. 

The model keeps matrix size and the number of grid 
elements small to keep the time necessary for solution 
within reason. The model concentrates grid elements in 
areas where the rate of change is rapid, or where the values 
are most important. Those areas are along the wing 
boundaries. An excellent method for concentrating the 
points, while keeping the total number of points at a 
minimum, is through cosine spacing. Piyice §7-5 > contains a 
sample of the layout. Notice that the rectangular elements 
do not model the planform exactly, but in the limit the 


difference will be negligible. 


Remote Velocity 


i 





Grid Model with Cosine Spacing 
Figure /-3 
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3. Matrices Used in VORTEX Program 
The program uses two identical ™ matric cee 
influence coefficients, [A] and [AA]. [A] is used to Wsomae 
for the ATI of the flat wing, and [AA] is used to solves 


the wing shape of the elliptically loaded wing.? 


—') dz 
al [#1 = fos} etrwian 





dz 
AA | jar} = | ” | (Elliptically Loaded Wing) 
x 


The SECTN matrix contains values for the wing 


sections, (chord, and width) as well as some of the final 
Eesules- Each row of the SECTN matrix corresponds tome 
wing section in the semi-span. The Output Variables 


section (H.l.a) describes the columns. 

The WING matrix contains wing coordinates required 
by the program as well as some of the final results. Each 
row of the matrix corresponds to an element in the wing 
semi-span. The Output Variables section (H.1.b) describes 


the columns. 


3The method used to develop wing shape will be 
discussed in) section, V Ui vs 
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4. Solution of the Set of Equations 


The subroutine GAUSS solves the N equations in N 
unknowns. The subroutine destroys [A] in the process but 
returns the solution vector (AI) in place of the right hand 
sides. 

5. Finding the Forces on the Elements 
The Kutta-Joukowski theorem states that the force 


per unit span acting on an element is. 


74 
Ak --___> poe 7 
—— = oey .. x oF ; 
Ay 


In this equation, Vere is the local effective velocity at 
the center of the element. Verr is defined in equation 7.8 
and Figure /-4. p is the density and AIF is the incremental 
Circulation around the element. iniserecreulation /CAL) is 
the same AIT contained in equation 7.1. AI can be termed the 
incremental circulation that occurs over the element. 

PirowmMOowemrnis Section, certain approximations and 
dimensional simplifications will be made. The first is the 
small angle approximation, where the sine is approximately 
Baual to the angle in radians, and the cosine is 
approximately equal to 1.0. | This approximation requires 
that angle of attack for the flat wing be small, generally 
less than 10 degrees. The approximation also requires that 
any wing slope on the elliptically loaded wing be small. 


This requirement is met by restricting the desired lift 
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coefficient to values Tess than apoumemue In the process 
of determining coefficients Of a Lili ts drag and moment, 
certain dimensional reference quantities arise. These 
quantities are density (pp), remote velocity (Vj), planform 
area (S), and average chord Cope ia performing dimensional 
analysis an arbitrarily unit value can be assigned to these 
berms. 
a. Finding ACy 

The distribution of AC, over the wing is 
desired, so the force on each element must be converted to a 
dimensionless pressure difference coefficient. Vote ee 


defined as. 


Vera 


Components of Effective Velocity (VaFrFP) 
Figuse /-4 
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Memis also sienificant to note that AY can be written as, 


A NP F 7.9 


since the bound portion of the horse shoe vortex is aligned 
with the y axis. Themresultesot the force cross product, 


memaction /./, 18s. 


AF > > 
— =p|(w—V_sina JAr: + V cosadlrk 710 
Ay [@ 6) wm 
Mmaecaneg p = 1, Ve. = 1, and incorporating the small angle 
approximation. 
= 
AF > + Tee 
—— = w—0)ari +ATR 
Ay 
This force is related to ACh. ACy Moe a SCalar. Father Chan 


a vector quantity, so the program uses the force component 


normal to the wing to compute ACy. iieae  eoupomene 1S. 
AF = AyAr yee 


Mepeime the f£oree non-dimensional gives a pressure 


difference coefficient, ACy. 
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Finally, after setting p» — | and=¥y =) AC, for an element 


es 


ry (adh Zee 
Neon = 7.14 
P Axdy 188 § 





This provides ste ACy for each element, which can be plotted 
versus chord, for each span-wise section. 
b. Finding the Lift, Drag and Moment 
Lift, drag and moment are also found from the 
vector force on each element. Equation 7.10 showed that the 
force on an element has components normal and tangent to the 
wing. Once more setting p = 1, and Vag = 1, the vecuen 


force on an element is. 
— ; > 
aF = |(w- sina JAP Ay? + cosadPay 7 alg 


These forces are oriented relative to the wing coordinate 


system, where Fi is parallel to the wing and k is 
perpendicular to the wing. Lift and drag are normal and 
tangent, respectively, to the remote velocity, ce The 


wing is at an angle of attack («) to that remote velocity. 


Using this information, the program transforms the forces to 
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wing tip and trailing edge (the nth elements) equal to a 
magerionwotuthe value of the [n-l]th elements. The program 


Pimiporees the condition by fitting a polynomial through the 


two points. See Figure 5-2. 


PeMote VEEREIry 
—_— —_———. 


N and N-]th Elements 
Figure 5-2 


4. Grid Spacing 
The program uses a non-uniform grid spacing, based 


on a cosine function. The elements near the wing tip are 


small. This spacing has the advantages noted earlier. 


C. PROBLEMS ENCOUNTERED 


In matrix form, the model has the form. 


« |lacn| = [2] - 





Mummies the influence coefficient matrix, and z° is the wing 


slope. The model specifies some values of AC), on the Lebe 


aS 


side of the equality and some values of wing slope (2°) on 
the right side of the equality. This complexity prevents 
use of the normal matrix solvers. As a result, the solution 
requires a complex matrix manipulation which is not readily 


understandable, or desirable for a teaching tool. 


D. CONCLUSIONS 

The evaluation of the integral over the element has 
definite advantages for coincident field and control points. 
The evaluation eliminates a strong singularity and should 
give a result using finite elements that is very close to 
the actual property. 

The Kutta condition and method used to enforce 


satisfaction is not optimal. Further studies are needed. 
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B. INTRODUCTION 

PGbRy Chiisel. HOLMewrote this program, as a tool for use 
ma AEZO35. Mies proccm@am is partial satisfaction of the 
requirements for the degree Aeronautical Engineer from the 
Naval Postgraduate School. The objective of the program is 
to provide a simple computer simulation of flow over a thin 
Wing with low aspect ratio. High aspect ratio wings are 
better served by a lifting line model. Many advanced 


programs using elaborate methods to model viscous effects, 


compressibility, boundary layer growth, and shocks, are 
available. None of them gives a good introduction to the 
field of computational aerodynamics. This program should 


fill that need at the Naval Postgraduate School. 

If you wish to skip the theory behind the program, go to 
section I for instruction on running the progran. 

There are a number of ways to model the source of 


Meeeces acting on a body surrounded by fluid in motion. 


fese inelude potential functions, vortex distributions, 
Mmeeulation distributions, and pressure differential 
cestributions. Each model is related to the other, and each 
has advantages and disadvantages. This program uses a set 


of horse shoe vortices to approximate the flow over a wing. 


SY 


Ga SYSTEM REQUIREMENTS 


System type IBM PC or compacreme 
Memory 256K, minimum 
Math Co-processor (8087) Recommended, Nn Ome 


required: 


Graphics card Recommen deja not 
required. 

Graphics program Xe Y graphing program 
like GRAPH ERS is 
recommended, n OMG 


required. 


D. CONCEPT DESCRIPTION 
l. The Horse Shoe Vortex 
This program, which models the flow by a series of 
horse shoe vortices, distributed over the wing, is an 
adaptation of the VORLAT program by Moran [Ret aiie@ He 
develops a program to find the strengths of a series of 
horse shoe vortices associated with a flat rectangular wing. 
The VORLAT program has its foundation in two-dimensional 
airfoil theory, and Moran [Ref.1] adapted the theory to 
wings of finite aspect ratio. 
The Users Manual will present a short description of 
the VORTEX program. For coverage of the VORLAT program, the 


foundation of the VORTEX program, consult Moran [Ref.1]. 
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a. Down-wash and the Relationship to Flow Tangency 
The velocity induced at a point in the xy plane 
by a horse shoe vortex, also located in the xy plane, can be 


determined from: 


e i i | V(x-x) +y-y,)° Sst 


= 1+ —— 
w (x,y) 4nly—y) (x—x) 4n(y—y,) 


(70) 
a 


me (xX = Xo), 


lh ; 
wen = 2 (| tie /4 
a b 
Gumeoee che velocity induced at (x,y), called a control point. 
Sets the strength of the horse shoe vortex, and (x,,Y,) and 
[/ee,ybh) are the corners of the horse shoe vortex. The 
center of the horse shoe vortex is a field point. 
TicmvettGCctiy at \a° point on a flat wing, Which 
Meeeat an angle of attack a, has components V,cos(a) and 


Veosin(a) as indicated in Figure 7-1. 


Veo 
Veo sinc 


Vao COSY FLAT WING 


Velocity Components on a Flat Wing 
Pogure 721 


ou 


In order to ensure the flow is tangent to the 
wing, the induced velocity or down-washt at that point an 
the wing must alse be @=Vosin (aye Each vortex on the wing 
contributes to the down-wash* at every control point on the 
wing. So, an equation for each control point, as a fume 
of all the field points@can bewyvreecen. An example of the 
equation for the control point (Xp +Vp) ise 


V sina~ s 
sina » See 0 Vix 3 
y 


When the control points are combined, the result is a Siete 
N equations in N unknowns. 
b. Placement of the Horse Shoe Vortex 
The bound portion of the vortex (that poreiem 
perpendicular to the onset flow) is placed at the 1/4 chord 
point of each grid element on the wing. Flow tangency 1s 


evaluated at the 3/4 chord point of each grid elemence 


lnown-wash is considered positive when its directions. 
downward. 


The contribution will be positive when the induced 


velocity is downward, or negative when the induced velocity 
ls upward. 
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ec. Satisfaction of the Kutta Condition 
Memesattisty s the, Kutta condition, the vortex 
strength must go to zero at the trailing edge. For wings 
and airfoils of zero thickness, flow tangency enforced at 
the trailing edge appears to be a corollary of the Kutta 
Bondition. 
2. Developing the Wing 
a. Planform Geometry 
It is necessary to develop wing geometry from 
aspect ratio, taper ratio and sweep angle. Aspect ratio is 
the span divided by the average chord. Taper ratio is the 
Grape chord divided by the root chord, and sweep angle is the 
mere between the leading edge and the y axis. Moran 
[Ref.1] uses a fixed root chord and varies the span as 
necessary to achieve the required aspect ratio and taper 


Bacio. The VORTEX program uses this approach, which works 


well. 
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Remote Velocity b/e 


° | 


Planform and Variables used in VORTEX Program 
Figwrem, -2 


Aspect ratio;  GAReeeess (b2/S). Span (b) is @eme 
distance from wing tip to wing tip, and area (S) IS.eime 
total planform area of the wing. 

Taper ratio, A, is the ratio of tip chord toni 


chord. 


tip chord © 7.4 


root chord Cc 
q 


Sweep angle delta, A, is the angle that the leading 


edge makes with a line drawn perpendicular to the "roeE 


Ciro ide 
b. Grid Geometry 
Using aspect ratio, taper ratio, and ~sweem 
angle, the subroutine SET85 develops a grid. The program 


stores the coordinates of all necessary points in the 
matrices WING and SECTN. Results of the program are also 
stored in WING and SECTN which are written to data files 


when the program finishes. 
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In a manner analogous to that used by Hough [Ref.3], 
Be G@ipexoGcices sane inset to improve accuracy of the 
results. The inset distance is one element wide. 

The model keeps matrix size and the number of grid 
elements small to keep the time necessary for. solution 
within reason. The model concentrates grid elements in 
areas where the rate of change is rapid, or where the values 
mee most important. Those areas are along the wing 
boundaries. An excellent method for concentrating the 
points, while keeping the total number of points at a 
minimum, is through cosine spacing. Figure 7-3 contains a 
sample of the layout. Notice that the rectangular elements 
Gomenot model the planform exactly, but in the limit the 


difference will be negligible. 


Remote Velocity 


—, ae 





Grid Model with Cosine Spacing 
Figure /-3 
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oe Matrices Used in VORTEX (oxoeram 


The program uses two identical matrices of 
influence coefficients, [A] and [AA]. {A] is used to solve 
for the Ar of the flat wing, and [AA] is used to sSolvVemieaa 


the wing shape of the elliptically loaded wing. ? 


lal” [] = [or] (Flat Wing) ie 





dz 
AA | jar} = | ” | (Elliptically loaded Wing) 
ie 


The SECTN matrix contains values for the wing 


sections, (chord, and width) as well as some of the £inal 
results. Each row of the SECTN matrix corresponds to a 
wing section in the semi-span. The Output Variables 


section (H.l.a) describes the columns. 

The WING matrix contains wing coordinates required 
by the program as well as some of the final results. Each 
row of the matrix corresponds to an element in the wing 
semi-span. The Output Variables section (H.1.b) describes 


the columns. 


3The method used to develop wing shape will be 
discussed in section VII.E. 


6 6 


& 


4, Solution of the Set of Equations 


The subroutine GAUSS solves the N equations in N 
unknowns. The subroutine destroys [A] in the process but 
jmeturns the solution vector (AI) in place of the right hand 
sides. 

5. Finding the Forces on the Elements 
The Kutta-Joukowski theorem states that the force 


per unit span acting on an element is. 


In this equation, eee is the local effective velocity at 
the center of the element. eee is defined in equation 7.8 
and Figure 7-4. p is the density and AI is the incremental 
circulation around the element. This een weculation (Ar): is 
the same AI contained in equation 7.1. AT can be termed the 
incremental circulation that occurs over the element. 

Throughout this section, certain approximations and 
dimensional simplifications will be made. The first is the 
small angle approximation, where the sine is approximately 
eoratl to the angle in radians, and the cosine is 
approximately equal to 1.0. its apprEomumation requires 
Smee angle of attack for the flat wing be small, generally 
less than 10 degrees. The approximation also requires that 
any wing slope on the elliptically loaded wing be small. 


This requirement is met by restricting the desired lift 
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coefficient to values less Gham apes In the procege 
of determining coefficients of lift, drag and momeme 
certain dimensional reference quantities arise. These 
quantities are density (p), remote velocity (Vs), planiueum 
area (S), and average chord (ey) in performing dimensional 
analysis an arbitrarily unit value can be assigned to these 
terms. 
a. Finding ACy 

The. distri bueion wae. ACy) over the wing is 

desired, so the force on each element must be converted to a 


dimensionless pressure difference coefficient. V off ee 


defined as. 


V0 W 


Vet 


Components of Effective Velocity (View 
Figure 7-4 
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Melis also Significant to note that AI can be written as, 


Ar=APry 7.9 


since the bound portion of the horse shoe vortex is aligned 
with the y axis. iiemmrestlt sot sthe —foree cross product, 


Paation 7./, 1S. 


AF ie 5 
fe ate w—V_ sina jAri + V cosadrk 7 1 
y 
pemeting p = 1, Veo = 1, and incorporating the small angle 
approximation. 
— 
AF > > Fait 
= 5 (w-ajar? + ark 
Ay 
Wes LTorce is related to AC,. AC is a scalar, rather than 


P P 


a vector quantity, so the program uses the force component 


normal to the wing to compute AC That component is. 


Dp: 
AF = AyAr tel 2 


Making the force non-dimensional gives a pressure 


difference coefficient, AC 


Pp ° 
AF 
AC = Viigo 
P 1 2 
ae Vi AxAy 


2S, 


Finally, after setting p=") and VZe-arm ACh for an element ~ 


ioe. 


ya Zar 
Nog =e 7 eee 
Po Oey Ax 





Tiis) provides. the ACy for each element, which can be plotted 
versus chord, for each span-wise section. 
b. Finding the Lift, Drag and Moment 
Lift, drag and moment are also found from the 
vector force on each element. Equation 7.10 showed that the 
force on an element has components normal and tangent to the 
wing. Once more setting p = 1, and Vom 1, the Vvecuan 


force on an element is. 
— ° > 
AF =|(w- sina JAP Ay? + cosaday k eS) 


These forces are oriented relative to the wing coordinate 


system, where i is parallel to the wing and Ke 
perpendicular to the wing. Lift and drag are normal and 
tangent, respectively, to the remote velocity (Vo). The 


wing is at an angle of attack («) to that remote velocity. 


Using this information, the program transforms the forcesiease 
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a coordinate system oriented to the remote velocity. For an 


m@dividwual element the lift and drag forces are; 


AL = Elemental Lift Force = (cosa) (AyAl) (cosa) — (w — sina) (AyATP) (sina) 


AD = Elemental Drag Force = (cosa) (AyAP) (sina) + (w — sina) (AyATP) (cosa) 


using a small angle approximation, 7 SA56 


Al. = Elemental Lift Force = AyAT — waAyAl = AyAT 


AD = Elemental Drag Force = AyAla + (w—a)AyAl’ = wAyAL 


Mimeme total lift and total drag acting on the wing is a 
Summation of the elemental lifts and drags. Using wing 


symmetry; 


einai =e AL 


semi — span 


Deal Dras = 28 AD 
semit—span 


The moment about the y axis for an individual 


eo 


element, with the accepted sign convention, is; 


AM = Elemental Moment = —(normal wing force) (x position of center of element) Td 


AM = Elemental Moment = (—AyAT)(x__) 


x 
cen 
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Total moment about the y axis is a SUmmatioOmees 


the elemental moments. Once again, using wing symmetry. 


M=TolalMoment=2 > AM 


semi—span 


62 Finding the Coefficients 


Non-dimensionalizing the total lift and drag gives, 








a g 
a. v°S 
g Pra 
7 aes 
D 
~— v°S 
g Pie 
where S = span x average chord = (b x c), p = 1 and VQ =e 


Non-dimensionalizing the total moment about the y 


axis gives, 


M 
C.,=-—-_ 
Mo ] sie Vin Il 
—OV oS 
2 fo 6] 
where c is the average chord. 
The -aerodynamic center X,, ass iiipercanme It is the 


chord-wise position about which the flat wing has zero 


moment. The program finds X27 enome 
xy _ ~ bletmoment 7.92 
ac total lift 
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Mais assumes that total moment is taken about the y axis, or 
x=0. tie jrornan = tainds the moment on the elliptically 
Ioaded wing, and uses the X., from the flat wing as the 


reference axis to find Cyac. 


E. TECHNIQUE FOR SPECIFIED LOADING 
1. Finding Camber and Twist from Loading 

Elliptic loading in a span-wise direction results in 
Minimum induced drag, and loading is proportional to 
Merewlation, At sufficiently high aspect ratio, a flat wing 
with elliptic area distribution generates elliptic loading 
but weighs more than a straight wing. Wineseewilth ‘straleht 
leading and trailing edges also cost less to manufacture 
than elliptic wings. So, a span-wise elliptic lift 
distribution is generated by slightly twisting the wing. 
The direct relationship between wing shape (camber. and 
Gwast) and circulation distribution makes this generation 


possible. 
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camber 





Camber and Twist on the Elliptically Loaded Wing 
Figure 7-5 


For a thin wing, elliptic load distribution "iia 
chord-wise direction requires an approximately parabolic 
camber shape. The program specifies an elliptic chord-wige 
load distribution and then finds the associated shape. 

2. Specifying the Lift Coefficient 

Specifying the form and amplitude of the ideal iigee 
distribution over the wing determines the form and amplitude 
of the wing shape. A load distribution of elliptic form is 
imposed in this case and the corresponding amplitude is 
fixed by the desired ideal wing lift soe fflerene Clas 

3. Maximum Desired Lift Coefficient 

The small angle approximation breaks down when large 
Phi) Coenriieiemts are specified: To prevent tims 
occurrence, the maximum desired lift coefficient is 


bBestricced, to -U.> . 
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ie Achieving the Desired Eete eCocttimcient 


The program scales a reference distribution of AC), 
which is elliptic in form, to achieve the desired total lift 
coefficient for the wing. 

wee Porces Moments and Coefficients 

Once the ACy (ict omenonew els Scaled for the 
elliptically loaded wing, the program finds forces, moments 
and coefficients in the same manner as for the flat wing. 
The only difference is in the coordinate transformation used 
to convert forces on the wing to lift and drag. The 
individual elements on the elliptically loaded wing are no 
longer at a uniform angle to the remote velocity. (bleak: 
change in angle must be taken into account when performing 
Siem coOrdinate transformation which gives lift and drag 
forces on individual elements. 


The program finds moment coefficient about’ the 


Gemoaynamic center, Cywac, from, 


where Xz. is found from the flat wing, c is the average 
enor d . Gworeisssenhe pitching moment of the elliptically 
iereed Wing about the y axis, and C;; is the (lift 


Seeeateient of the elliptically loaded wing. 


J 


Be FLOW CHART FOR MODULES 


SET WING SET GRAPHIC VIEW GRAPHIC END 
PLANF ORM DISPLAY DISPeAr PROGRAM 
COMPUTE 





CORREICIENTS 
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INPUT VARIABLES 


There are eight input variables. They are: 
Aspect Ratio. AR, defined as b2/S, where b is the 
total span and S is the total area, can have any 
positive value less than 20.0. Higher aspect ratios 


should use lifting line theory. 

Taper Ratio. LAMDA, defined as C;/C,, where Cy is the 
wing tip chord and Cy, is the wing root chord, can have 
any positive value between O and l. 

Leading Edge Sweep. DELTA, defined as the angle in 
degrees that the leading edge of the wing makes with a 
line perpendicular to the remote velocity, can have 
any positive value between 0 and 60. 

Angle of Attack. ALPHA, defined as the angle in 
degrees made between the flat wing and the remote 
velocity, can have any positive value between O and 
12O®. 

Number of elements in a section. NX, defined as the 
integer number of elements in a chord of the wing, can 
have any value between 1 and 10. 

Number of elements in a semi-span. NY, defined as the 
integer number of elements in a semi-span of the wing, 
can have any value between 1 and 10. 

Desired lift coefficient. CLDSRD, the desired lift 
ecetficient for the elliptically loaded wing, can have 


any value between O and 0.5. 


Dd 


H. OUTPUT VARIABLES 


a le Tabular Date 


When complete, the program writes four tabular data 


fives: One file, the SECTION.MAT file, ¢eoentaine waa 
each section in the wing semi-span. Another, the WING 
file, contains data for the individual grid elementsmiam 


wing semi-span. The third, the FLAT.DAT file) Sehere 


coefficients for the flat wing. The fourth, the CAMBER 


£2 de. shows the coefficients fon the cambered 
elliptically loaded wing. 


a. SECTION .MAT 


Lou 


. Mad 


the 


Each row, or line, in the file represéma@eue 


section of the wing semi-span. Each column represents a 


different variable within the sec e1one 


ye 
Remote Velocity 


ne 


chord 


A Single section in the Wing Semi-Span 
Figure /-6 
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mhe columns are: 


ie Section chordm@leng th 
temeoection lift scoefficient per unit span for the flat 
wing. The integral of these values from tip to tip is 


the total lift coefficient. 

wee cection drag coefficient per unit span for the flat 
wing. The integral of these values from tip to tip is 
the total drag coefficient. 

4. Section moment coefficient per unit span about the 
aerodynamic center for the flat wing. The integral 
of these values from tip to tip is the total moment 
coefficient about the aerodynamic center. 

Be): ; section Xgc, aerodynamic center relative to x=0, for 


the flat wing. 


SO. Section delta y 
Mumeoecetron. Litt coefficient per unit span for the 
cambered wing. The integral of these values from tip 


Mmm ic the total tiftt coefficient. 

8. Section drag coefficient per unit span for the 
cambered wing. The integral of these values from tip 
Pomcipelswthe total drac coefficient. 

Eee oection moment coefficient per unit span about the 
aerodynamic center for the cambered, or elliptically 
loaded wing. The integral of these values from tip 
to tip is the total moment coefficient about’ the 


aerodynamic center. 


1 


10. Section twist in degrees of the cambered wing. ~ 
Positive angle of twist is nose up. 

The program arranges the rows from root to tip, with the 

root section on the first line and the tip section One 
last line. 

be WING.MAT 

Each row, or line in the file, representsiiaa 

individual element in the wing grid. Each column represents 

a different variable for that element. Figure 7-7 shows the 


location of some variables on a representative grid element. 


a a es 
R Veloci 
suse eens JAW 


é] 


a 
(Xvi? 


(XB, Yep (POree® 


; (XP, YP) (shape? 


Cx, ya 


— delta x 


A Single Element in the Wing Semi-Span 
Figure 7-7 
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The columns are: 


ke 


nO. 


He... 


rE? . 


ee 


14. 


Joa 


EG 


li 


Eo . 


Ae? 


ScquecmtraimmEerser, “1adentitiler of the element, also 
called IJ in the progran. 

Integer y position of the element. iMais@wae Ehe wine 
BOOtCemNiels at the tip. 

Integer x position of the element. eas at. the 
leading edge, NX is at the trailing edge. 

x position of center of element. 

Yepesition of eenter of element. 

MEPociWwtclolew Or Worse. Shoe vortex, a field point. 

Blank, or zero. 

ya position of one corner of horse shoe vortex. 
VOepOStebom= er One corner of horse shoe vortex. 
Ppecoaiiline edge of the element, a control point used 
in computing wing shape or AT. 

Mc cllecimmor "the clement, a control pointm used in 
computing the forces acting on the wing. 
PPmceilcerwOLuwtmeuelement, a control point. 

Delta x, chord-wise dimension of individual element. 
Delta y; span-wise dimension of individual element. 
MWetcwine ine remental circulation (AI). 

Eilat wing ACy. 
tampered awite “ineremental cireulation (AIT), adjusted 
Ome GT ret - 


Cambered wing ACy 


Cambered wing slope (dz/@x). 


onl 


ZO. 


pag ier 


LP ee 


2a 


Some 


each 


Cambered wing height (z). 

Cambered wing incremental circulation (Al), if )G@igee 
were = 1.0. 

Down-wash (w) at center of flat wing element due to 
trailing filaments. 

Down-wash (w) at center of cambered wing element due 


to traidame £1 Vamentee 


columns are identical. Separate columns were used for 


variable during the development phase to promote ease 


in modification. The rows, or lines, of the WING.MAT file 


are arranged as follows. 


The 


Remote Velocity 





Grid Element Numbering on the Semi-Span 
Figure 7-8 


first row is the leading edge element in the root 


section. The next element is immediately behind the leading 


edge and so on. The last row is the trailing edge element in 


the wing-tip section. 
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Mes 


eee ELATSDAT 


VimoGeetatleommcontarins coetficients and results 


mae flat wing. 


FLAT WING 


Gr 

CD 
eD/CE2 
CMAC 
XAC 


ALPHA 
CLDSRD 
Oey Ale 


8 
8 


wis 


S 


An example is shown. 


SAAT 
POUGo iG 
PeUGo i 

. 000000 
2799285 
Ig 
. 900000 
Ne 


330000 


000000 


. 000000 
. 200000 


CLa, Lift Curve Slope 
per/Deg= 


d. CAMBER. DAT 


.035544 


MMe imMieNcrMCGiicatmSueGOect ficients and results 


the cambered, 


shown. 


CAMBERED 


GL 

CD 

en /GL2 
CMAC 
AR 
LAMDA 
DELTA 
NX 

NY 
ALPHA 
CLDSRD 
ELLIP 


WING 


YES 


8 3 


Out ert ientrcally, Loaded, 


poo cam 
009040 
L215) 
0G 27 32 
SECO, 
. 900000 
000000 


000000 
. 200000 


wing. 


An example 


for 


£ o'r 


is 


aie SAMPLE PROBLEM 


A sample problem will illustrate use of the |VORsMas 
program. A wing planform with the following characterise 


will be analyzed. 


ASpect  racro fighes eh 
Taper ratio Or. 
Leading Edge Sweep 25 degrees 
Angle of Attack > degrees 
Number of span elements 8 


Number of chord elements 8 
Desired Lift Coefficient 90.2 


The planform looks like Figure 7-9. 


Ke 0.5000 —=y 


ao | 


ReMOtes Jelegiiry 0.5000 


Planform used in the Sample Problem 
Figures 22 
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Ik. Starting the Program 


After turning the computer on go to the DOS prompt, 
which generally looks like this. 
oe 
Change the program to the VORTEX directory by typing 
CD\VORTEX [return] 
the screen should now look like. 
ev ORTEX : > 
To start the program, type 
VORTEX [return] 
The program will start and the following menu will be 


displayed. 


MAIN PROGRAM MENU 


Set Wing Planform, compute Coefficients and Shape 





Set Graphic Display 


Use this option to set up 


View Graphic Display VOourT wine Dplkantorm. It 
will also compute al! the 
End the Program coefficients you need, 


Lift, [rag, Moment, as 

well as the share, if you 
eee EE want an ellaptic load 

Gistrinition. 





Matinee rogram Menu 
Figure 7-10 


2. Main Program Menu 
There are four options in the MAIN PROGRAM MENU. 
All coefficients and results are computed in option 1. 
Whenever the planform is changed, option 1 must be selected 
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to recompute the coefficients. The coefficients and 1eswiaaes 
are saved to data files after being computed in opt tome 
The data files allow you to use options 2 and 3) as (Geet 
desired without re-computing the results for the plantomm 

Help menus are available. Press Fl to topple wae 
help menus on or off. 

3. Planform and Coefficient Menu 

The program saves planform variables and results 

from the most recent solution and shows them in the PLANFORM 


and COEFFICIENT MENU. The menu looks like Figure 7-11. 


Re mt ee ee 


PLANFORM AND COEFFICIENT MENU 









1; ASpect Ratio .. 224444454464 )5)22 ee 1.330000 
2s. Taper Ratio Wis c¢ sitet. 1 eee -500000 
3. Leading Edge Sweep Angle (degrees) ..... 25.000000 
4. Angele of Attack (decrees) S27. 5.000000 
S. Number of Elements ina Secretion ee 8B 

6. Number of Elements in a Semi- Span ee 8 

7. Compute Wing Shape with s boads. ae YES 

8. [Lesired Lift Coetiicienk 4... ee . 200000 
9. Compute Cofticients £or this sl lannomun 
10. Return to Main Menn 





Aspect ratio has values 


between O and 10. 





Planform and Coefficient Menu 
Raber ee 7 le 


Options 1, 2, and 3 change Aspect Ratio, lapes 
Ratio, and Leading edge sweep angle, respectively. 

Option 4 changes the angle of attack. 

Option 5 changes the number of chord-wise (jou 


elements, and 
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Option 6 changes the number of span-wise_ grid 
elements in a semi-span. 

Option 7, Compute Wing Shape with Load, is YES or 
NO. Because the time required to compute wing shape is 
extensive, that option may not be desirable for all 
problems. 

Option 8 changes the desired lift coefficient, Cy]j 


of the elliptically loaded wing. 


Option ) computes the CoemriLcilents . The 
eoetficients are written to rants files, SECTION.MAT, 
WING.MAT, FLAT.DAT, and CAMBER.DAT. Section. H of this 


thesis contains a description of the elements in the files. 


Option 10 returns to the Main Program Menu to view 


ne results. 
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4. Graphics Program Menu 
Choosing selection 2 of the Main Propramegtegm 
displays the Graphics Program Menu. An example is shown in 


Fleer eae? 






GRAPHICS PROGRAM MENU 






CAMBERED 
> 





Type Of Wing. 14:4. e nee 
Section Number 22.4... eee 
Return to Main Menu 







Spanwise Plots 






AGE Malay 

ACW ay, 

d(CMAC )/dy 
Camhered Wing Twis 






Chordwise Plot|} to the Main Program Menu, 
where you can view the 
Delta Cp graph you have set up 
Wing Shape here. 





ee ee 


Graphics Program Menu 
Figure 7-12 


Option l changes the type of wing that wile 
displayed, either flat or cambered. 

Option 2 changes the section of the wing which is 
plotted with options 8 and 9. 

After setting options 1 and 2, select the desamged 
plot from options 4 through 9. The necessary data files 
will be prepared and the program will return to the Main 
Program Menu, where the selected plot can be seen using the 
View Graphic Display option. Examples follow. 
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FLAT WING 
Gia = .1693, AOA = 5.00 deg 


d(CL)/dy 


Oo 
oh 
nN 


0.50 


0.24 


= 
OD 


0.06 


0.0 0.1 G2 O-2 0.4 8s, 


| SPAN—WISE COORDINATE Y | 


FLAT WING d(CL)/dy vs SPAN 
Figure 7-13 


| CAMBERED WING | 


CLi = .1993, CLi = a4 0, 


CAMBERED WING d(CL)/dy vs SPAN 
Figure /-14 


FLAT WING 
CD = .0062, AOA = 5.00 deg 


d(CD)/dy 


=9):020 


0.020 


0.012 


OO12 





SPAN-—WISE COORDINATE Y 


FLAT WING d(CD)/dy vs SPAN 
Figure 7-15 


| CAMBERED WING | 
CD = .0090, CLi = 220 


d(CD)/dy 


0.020 


0.016 


0.012 


0.008 


0.004 


0.000 


0.0 0.1 g.2 G3 0.4 0:5 : 
| SPAN—WISE COORDINATE Y | 


CAMBERED WING d(CD)/dy vs SPAN 
Figure 7-16 
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| FLAT WING | 
CMAC = .0000, AOA = 5.00 deg 


d(CMAC) /dy 


G.011 


0.0070 


0.0030 


! 
= 
Q 
o 
oS 


—0.0050 


—0.0090 


0 


oO 





| SPAN—WISE COORDINATE Y | 


FLAT WING d(CMAC)/dy vs SPAN 
Figure 7-17 


| CAMBERED WING | 


d(CMAC) /dy 


CMAC = —.0625, Cli = 20 


—0.000 i =) Sn on) et 
0. 0.1 O72 0:5 0.4 Oro 


SUG lS 
—0.056 
—0.054 


=—@.07 2 


| SPAN~—WISE COORDINATE Y | 


CAMBERED WING d(CMAC)/dy vs SPAN 
Figure 7-18 
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J. METHOD FOR CHANGING LOADING DISTRIBUTION 

This is an optional section and may be skipped without 
loss. Elliptic loading is generated in the  SET85 
subroutine. The elements which contain the specified 
loading are WING(IJ,21). If you wish to change the loading 
from elliptic, realize that the grid elements are not 
uniform in size. 

When generating a ACy function, consider the shape. The 
load at the wing tips and trailing edge must go to zero. 
The load at the leading edge should preferably be zero to 
evOld a Singularity at that point. 

If you change any of the subroutines, the program must 
be re-compiled with a FORTRAN compiler to include _ the 
Smeameges into the VLAT85.EXE file. The original version of 
the program was compiled by a MICRO-SOFT 4.2 compiler using 
the following commands. 
miytPc VLAT&85.FOR 
and 
Reyer Pc GRAP85.FOR 
If you have a math co-processor, the FPe option allows the 


program to use the co-processor but does not require one. 


K. METHOD FOR CHANGING THE NUMBER OF GRID ELEMENTS 

This is also an optional section and can be skipped. iB el 
Mess Original form, the program works with up to 100 grid 
elements in the wing semi-span. The program also has a 
maximum of 10 chord sections. Fewer grid elements or chord 
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sections will run without modification. If you desire more 
than 100 grid elements or 10 chord sections, you must modify 
the program. The changes are not extenSive, but require 


that the programs be recompiled, as described in section J 


above. 
Increase grid size in square increments. Modify the 
matrix elements as follows. N is the square dimension and 


corresponds to NX or NY. 

For example, if you want 12 wing sections or 12 eee 

elements in each chord section, make the grid 12x]12 ie 
elements). In this example, N = 12, and (N*¥N)+1 = 145. The 
matrices that must be changed are listed below. Tire 
required dimensions are shown. 
ACCN*¥N)+1, (N¥N)+4+1) 
AAC CN*¥N)+4+1, (N*¥N)4+1) 
AS (C(N*¥N)+4+1, (N*¥N)+41) 
WING( (N*¥N)4+1,23) 


SECTN(N+1,10) 


Oe 


These are the line numbers in the programs’ which 
eontain matrix variables that must be changed. 


Variable 


Program A() AAC) AS() WING() SECING } 
VLAT85 76 a, Vy 74 74 
LOAD85 56 ay a7 54 54 
SHAP85 54 3 a) DU 52 
GAUSS Ih, - - z a 
MULT1 ie, 20 AG 17 Jy) 
iw LT 2 IS 20 20 iy 17 
PuLT 3 Jo) Ze 746) Ly iy 
GRAP85 Dic 53 3 50 50 
DNWS85 - - - a a 
SET85 - - - 69 69 
Before increasing the grid density, cComsider the 
tradeoff of memory requirement and solution speed. The 10 X 
10 grid contains 100 elements and 100 unknowns. A SEZ x 12 
grid contains 144 elements and 144 unknowns. A 20 x 20 grid 
contains 400 elements and 400 unknowns. More than 4 times 


the number of computations are needed to solve a 400 x 400 
Mermsus a 100 x 100 matrix. You should also increase 
precision if you increase the number of grid elements. All 
elements will be smaller and computer accuracy may be less 


than the distances between points in the grid. 
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APPENDIX A. SOURCE CODE LISTINGS FOR VORTEX PROGRAM 


VLAT85 


tek kakkiek kok kkk kkk kkk ok ek seat ok ake ok kak ok ak ok ak ak ok ak ak ak ak ok ak ok ak ak ak ak ak ak ok ok ok ak ok ok sh ak de ak ok ok EK ak 


* * 
* PROGRAM VLAT85 * 
* * 
* THIS IS THE MAIN DRIVER FOR A NUMBER OF SUBROUTINES THAT ARE * 
* LINKED TOGETHER TO FORM A PROGRAM THAT COMPUTES THE FORCES * 
* ON A PLANFORM IN INVICID, INCOMPRESSIBLE FLOW. * 
* THE SUBROUTINES CALLED ARE LISTED HERE AND AT THE END OF THE * 
* PROGRAM AS INCLUDE FILES. * 
* * 
* THIS PROGRAM IS A HIGHLY MODIFIED VERSION OF A PROGRAM IN * 
* ‘AN INTRODUCTION TO THEORETICAL AND COMPUTATIONAL AERODYNAMICS’ * 
* BY JACK MORAN, WILEY AND SONS, NEW YORK, 1984, FOUND ON PAGE * 
* Soke * 
* * 
* PROGRAM VARIABLES ARE: * 
* AR : ASPECT RATIO = 
* LAMDA : TAPER RATIO * 
* ADELTA : LEADING EDGE SWEEP ANGLE, IN DEGREES * 
* DELTA :; LEADING EDGE SWEEP ANGLE, IN RADIANS * 
* NX : NUMBER OF GRID ELEMENTS IN A SECTION * 
* NY : NUMBER OF GRID ELEMENTS ALONG A SEMI-SPAN * 
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+ 


ALPHA 
CLDSRD 


EEL Ie: 


NS 27Gs) 


WING( ): 


: ANGLE OF ATTACK, IN DEGREES 


DESIRED LIFT COEFFICIENT WITH ELLIPTIC Re wis 
CHARACTER FLAG TO COMPUTE THE SHAPE FROM AN 
ELLIPTIC LOADING 

Ss. Pel 


NX*NY, THE TOTAL NUMBER OF EQUATIONS 


: A COEFFICIENT MATRIX, USED FOR THE FLOW TANCENG 


GETS TRANSFORMED IN GAUSS 


THE SAME AS A BUT NOT TRANSFORMED IN THE PROGRAM 


: ANOTHER COEFFICIENT MATRIX, USED FOR THE DOWN-WASH 


AT THE CENTER OF THE GRID ELEMENT. REQUIRED TO 

GET THE FORCES ACTING ON EACH GRID ELEMENT. 

A LARGE ARRAY HOLDING VARIOUS IMPORTANT COORDINATES 
FOR EACH GRID ELEMENT, AS WELL AS THE RESULIC@He 
EACH GRID ELEMENT. A COMPLETE DESCRIPTION OF EACH 
ELEMENT IN THE ARRAY IS CONTAINED IN THE SUBROUTINE 


SET85 


SECTN(): A LARGE ARRAY HOLDING SOME DIMENSIONS FOR THE 


Lin 


LOUT 


JOUT 


KOUT 


LOUL 


SECTIONS, AS WELL AS THE COEFFICIENTS FORVEAGH 
SECTION. A COMPLETE DESCRIPTION OF EACH EEE 
IN THE ARRAY IS CONTAINED IN THE SUBROUTINE SET85. 
SET TO 5, FOR KEYBOARD INPUT 

SET TO 6 FOR SCREEN OUTPUT 

SET TO 10 FOR WING ARRAY INPUT AND OUTPUT 

SET TO 11 FOR SECTN ARRAY INPUT AND OUTPUT 


SET TO 12 FOR LAST PARAMETER INPUT AND OUTPUT 


104 


* MOUT 

* NOUT 

* IMOD 

* 

* 

* ICHOIC 
* ESC 

* ITYPE 

* ISECT 

* 

* SUBROUTINES 
* DATS . 
* CDIT85 , 
* LOAD85 , 
* 

* SHAP85 , 
* 

* DELAY 


SET TO 13 FOR FLAT PLATE WING COEFFICIENTS 


SET TO 14 FOR CAMBERED WING COEFFICIENTS 


: A FLAG. IF THE PARAMETERS ARE UNCHANGED FROM THE 


MOST RECENT CALCULATION OF COEFFICIENTS, SET = 0 


IF ANY PARAMETERS ARE CHANGED, SET = 1 


: MENU INPUT VARIABLE 


b] 


THE ESCAPE CHARACTER (27) 
WNPEGER. 1 IF FLAT, 2 IF ELLIPTIC 


INTEGERS SECTION NUMBER TO BE PLOTTED. 


CALLED ARE: 

FORTRAN, GENERAL DATA ENTRY, REAL 

FORTRAN, GENERAL DATA ENTRY, INTEGER 

FORTRAN, COMPUTE LOAD DISTRIBUTION FROM A GIVEN 
WING SHAPE. 

FORTRAN, COMPUTE SHAPE FROM A GIVEN LOAD 
DISTRIBUTION. 

FORTRAN, SLOW DOWN THE PROGRAM SO ERROR MESSAGES 


CAN BE SEEN BEFORE BEING ERASED. 


* 


* 


HHH HHH KKK KK KH KKK KKH K HK KKK AK KEK HK KEK EK HE KEK KEK KEK KEKE KEKE KEK EEK EREREEREENK 


PROGRAM VLATS85 


REAL LAMDA 


CHARACTER*3 ELLIP 


CHARACTER*1 ESC 


COMMON PI,B,WING(101,23), SECTN(11,10) ,AR, LAMDA, DELTA, NX,NY, 


OS 


* ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101), NEQNS 
COMMON /CAM/AA(101,101), AS(101,101) 


COMMON/FILS/TIN, IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 


IMOD = 1 
AR = 220 
LAMDA = 1.0 


ADELTA= 10.0 


NX = 4, 
NY = d, 
ALPHA = 5.0 


eo 
SS 
S 
ke 
rg 
| 
H 
tH 
Ww 


CLDSRD= 0.5 
[PYRE =v 
{SECT — 1 
XAC = 0.0 


PI = 2.0*ACOS(0.0) 
ESC = CHAR(27) 


KAKKHKH HHA HAH KH HHH HW Wk Wk Hw Wk ak Hk ake ak ok ak ok ak sk ok ak ok ak ak ok ak ak ok ak ok ok ak ak ok ak ok ok ak ak ok ak ok ok ak ak oh ok ak ak a KK 


* ¥ 
* SET INPUT/OUTPUT VARIABLES AND OPEN THE NECESSARY FILES * 
* MOUT AND NOUT ARE NOT OPENED UNTIL THE COEFFICIENTS HAVE BEEN * 
* COMPUEED: a 
* * 


Haka ak ak ak ak He ok ak ak ok ok ok ok ak ak ob ok ak ok ok ak ak ok ak ok ok ak ok ok ok ob ak ob ob ak ak ok ok ok ak ok ob ak ak ok ok ok ok ak ok ok ok ok ok ak oh ok hk KK EK KH 


IIN = 5 


LG 


Bed 


IOUT = 6 


JOUT = 10 
KOUT = 1l 
LOUT = 12 
MOUT = 13 
NOUT = 14 


OPEN(UNIT = IIN) 


OPEN(UNIT = IOUT) 
OPEN(UNIT = JOUT, FILE='WING.MAT’ ) 


OPEN (UNIT 


KOUT, FILE=’SECTION.MAT’ ) 


OPEN(UNIT = LOUT, FILE=' PARAMS.DAT’) 

10 FORMAT (16/3(F6.2,/),2(13/),F6.2/A/F6.2/12/12/F6.2) 

30 FORMAT (23F10.6) 

40 FORMAT (10F10.6) 

50 FORMAT (' THAT IS NOT A VALID OPTION. ENTER AN OPTION BETWEEN 1’ 
* ' AND 11‘) 

60 FORMAT (' ',A,A,F10.6) 

62 FORMAT (' ‘',A,A,13) 

64 FORMAT (' ',A,A,A) 


KHKEA KE AKA KEK KEK KEK KE KKK K EHH KKK KEK KHEKWK EH HEH HEH W HW WH WH WHE WW W ME ke 9k 9k ak Hk Hk 


* * 
* READ THE LAST SET OF PARAMETERS, IF UNCHANGED, READ WING AND * 
* SECTN DATA - 
7 * 


HHHHHAH KKK AK KH WK HHH HHH HHH HAH HHH HW HK HW HW Hk HH" KH Hk ve We ak Hk ok sk ok ok sk ok ok ak ok ak Hk ok ok ak ok ok ok ok ok IK ak 


READ (LOUT,10,ERR=70,END=70) IMOD,AR,LAMDA,ADELTA ,NX,NY,ALPHA, 
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* ELLIP,CLDSRD, ITYPE, ISECT, XAC 
70 DELTA = ADELTA*PI/180. 
B = AR*(1.0+LAMDA) /2.0 
NEQNS = NX*NY 
IF (IMOD.EQ.0) THEN 
READ (JOUT,30) ((WING(I,J),J=1,23), I=1,NEQNS) 
READ (KOUT,40) ((SECTN(I,J),J=1,10) ,I=1,NY) 
ELSE 
ENDIF 


Fe oe Xe oe ak ok ob Xe ak ok ok oe ak ok] ok ok ok ak ok ok ak ak ok ok oe ak ak ok ok ak ok ob ok ak ak ok ok ak ak ok oe ak 2k ok ok 2k He He oe IE ae ak 9] oe ae ak ae ok ok ae oe Oe aE EK 


* * 
* CLEAR THE SCREEN AND WRITE THE MAIN MENU. * 
* THIS USES THE ANSI.SYS CONTROL CODES TO CLEAR THE SCREEN AND * 
* POSITION THE CURSOR. IT ALSO USES FLASHUP WINDOWS FOR THE * 
* MENU «PROMPTS. ~ 
* * 


RIK KKK HHH KKK HH MH KH HHH HH HHH HHH HHHHMH KH HK HHH HH KH HH HH HH HH HHH HK NHK KW HK HH KH HW HK 
100 WRITE (lOUTs = Ese (29 

WRETE “(CLOUT | +) C=Ataey as 

WRITE (IOUT,*) ' W=LOAD/’ 

WRITE (IOUT,60) ESC,’(09;69H’ ,AR 

WRITE (IOUT,60) ESC,’[10;69H’ , LAMDA 

WRITE (IOUT,60) ESC,’[11;69H’ ,ADELTA 

WRITE (IOUT,60) ESC,’[12;69H’ , ALPHA 

WRITE (IOUT,62) ESC,’[13;69H’ ,NX 


WRELE S(LOUT262) a ESGaikL4 ; OSes 
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WRITE (IOUT,64) ESC,’[15;69H' ,ELLIP 
WRITE (IOUT,60) ESC,'[16;69H’ ,CLDSRD 
NEQNS = NX*NY 


KKK KKK KKH KK KKK KKK HHH KKK HK KKK KKK KEK KKK KK KKK KKK KB KK KEK KEKE KEK KE KEKE KKK ERE 


* * 
* READ THE OPTION. IF OUT OF RANGE, PRINT ERROR MSG, AND START * 
* OVER. * 
* * 


KEKKKHKKEKK KKK KE KK KH KHAH KK KKEK KKH HK HHH HHH HHH KKH HHH HHH HHHH HK HHH KKH HH KHK KKK KKK 
110 READ(IIN,*,END=110,ERR=110) ICHOIC 
WRITE (IOUT,*) ESC,'[2J' 
WRITE (IOUT,*) ESC,’ [00;00H' 
IF ((ICHOIC.LT.1).OR.(ICHOIC.GT.11)) THEN 
WRITE (IOUT, 50) 
CALL DELAY(6) 
GO TO 100 
ELSE 
ENDIF 


KHKKAKKEKEKK EK KEK KEK KE KK EKER KER EER EERE REE ERK EKER EKER EEK KEE EK KEKE EE LEK EK 


* * 
* MAKE ALL CHANGES TO PARAMETERS AS DESIRED. * 
* VALUES PASSED TO CDAT AND CDIT ARE LOW, HIGH, VARIABLE, AND * 
* NAME OF THE VARIABLE. IMOD IS RETURNED WITH VALUE 1. * 
* * 


KRKEKEKKE KKK KKK KEEN KKK KEKE KER EKER ER ERE RE RRR RRR REE RE RR ERE KER ERE REE ERE EK EEK 


oe ECHOIC EQ 1) CALL CDAT85(0.0,20.0,AR, 


Oo 


* ‘ASPECT RATIO ' IMOD) 
IF (ICHOIC.EQ.2) CALL CDAT85(0) 0) Ion LAM 
* ‘TAPER RATIO ' IMOD) 
IF (ICHOIC.EQ.3) THEN 
CALL CDAT85(0.0,60.0,ADELTA, 
* ‘LEADING EDGE SWEEP ANGLE IN DEGREES ‘,IMOD) 
DELTA = ADELTA*PI/180. 
ELSE 
ENDIF 
IF (ICHOIC.EQ.4) CALL CDAT85(0.0,10.0,ALPHA, 
* "ANGLE OF ATTACK IN DEGREES ' IMOD) 
IF (ICHOIC.EQ.5) CALL CDIT85(0,10,NX, 
* ‘NUMBER OF ELEMENTS IN A SECTION ' IMOD) 
IF (ICHOIC.EQ.6) CALL CDIT85(0,10,NY, 
* ‘NUMBER OF SECTIONS IN A SEMI-SPAN  ',IMOD) 
IF (ICHOIC.EQ.8) CALL CDAT85(0.0,0.5,CLDSRD, 
* ‘DESIRED LIFT COEFFICIENT ' IMOD) 
IF (ICHOIC.EQ.7) THEN 
IMOD = 1 


IF(ELLIP.EQ.'YES') THEN 


ELELP = NOs 
ELSE 
ELLIE T=. YEo 
ENDIF 
ELSE 
ENDIF 


hike, 


sea I II I Ra ae ae ok oe ak ok ak ok ak ok ok ak ok ok ak ok ak ok ak ok Xk akc ok ak ok ok ok ok ok oh ak ok ok ok 9k ake ok ok ok oh ok ak ok ok ok ok ok ob a ok OE oe Oe Oe ae Ok OK OE 


* 


* 


* 


* 


* 


RUN THE MAIN PART OF THE PROGRAM. 


COMPUTE THE LOAD, FROM THE 


FLAT PLATE SHAPE, AND THEN THE SHAPE, FROM THE ELLIPTIC LOAD 


ITF THE ELLIP FLAG IS SET. 


* 


* 


se akok ok ake ok ok ok ak ok ok ak ak ak ak ok ak ok ae ak ae ok ak ak ok ae Kok ae aE ak ae oe ak ok ok ak ok ok ak ak ok ak ok ake ak ok ok ok ak ok ok ok ok ok ok ok ak ok ok oe ak ok Ie Oe Oe Oe IE IE 


HS 


Pie LGHOTG; EQ. 9) THEN 
IMOD = 0 


WRITE (IOUT,*) ESC,'[2J' 


WROE, (TOUT. *) ESC, [00: 00H’ 


WR (TOU «*)-! C=ALL/’ 
WRITE (IOUT,*) ' W=WAIT/’ 


CALL LOAD85 


IF (ELLIP.EQ.'YES’) CALL SHAP85 


ELSE 


ENDIF 


FHI HIE FA EE BMW HW HB BW HWM HH WH "ke Me EN He Me Ok OE IE ak ok oe ae ok ok ok ok ok ok Hk ok ok ok ok ah ok ok ok Wh ak ok ok oe ok 


* 


* 


* 


* 


OPTION 10, RETURN TO THE MAIN MENU. DO A NORMAL END. WRITE THE 


PARAMETERS AND DATA TO FILES AND CLOSE THEM. 


ae 


¥* 


WHR KKK KK KK KE KKK KKK KK KKK KKK KKK KHER HH HW HW WW HWW WH WW FH HWW BE EE EE EEE EEE EE 


IF (ICHOIC.EQ.10) THEN 
REWIND (JOUT) 


REWIND (KOUT) 


ale 


REWIND (LOUT) 
WRITE(JOUT,30) ((WING(I,J) ,J=1,23),I=1,NEQNS) 
WRITE(KOUT,40) ((SECTN(1I,J) ,J=1,10) , I=1,NY) 
WRITE(LOUT,10) IMOD,AR,LAMDA,ADELTA,NX,NY,ALPHA, ELLIP, 
* CLDSRD , ITYPE, ISECT, XAC 
CHOSE (IOUT) 
CLOSE (KOUT) 
CLOSE (LOUT) 
GOTO 1000 
ELSE 
ENDIF 
GO" 10.4100 
1000 WRIEE (1lOUT + )2ESO {20 
WRITE (IOUT,*) ESC,’[00;00H’ 
WRITE (IOUT,*) ' C=ALL/' 
WRITE (IOUT,*) ' W=MAIN/’ 
END 


KKKKKEKKEKEKAA KEK EEK EK KEK KEE KEK KEE ERK EK KEKE KK EK KEE KEE EEE EEE KEKE KEE KEENE EK KENT ENK 


* * 
* THESE ARE ALL THE SUBROUTINES THAT ARE LINKED, * 
* * 


HK KH HH HA AAA HH HAHAHAHAHA KAKA A KKK KKK KKK KEK HHA KH KKK KHANH HHH HH 
SINCLUDE: 'SET85.FOR' 
S INCLUDE: 'DNWS85. FOR’ 
SINCLUDE: ‘GAUSS. FOR’ 


S INCLUDE: 'CDAT85. FOR’ 


Tale 


—; 


SINCLUDE: 
SINCLUDE: 
SINCLUDE: 
SINCLUDE: 
SINCLUDE: 
SINCLUDE: 


SINCLUDE: 


s6DITS >. hOR: 


‘LOAD85 . FOR‘ 


‘SHAP85. FOR‘ 


‘DELAY. FOR’ 


TET: FORe 


-MULE2Z. FOR. 


MUIET3 FOR) 


3 


LOAD85 


KKKKKEKKEKE KERR KEKE EER EKER ERREKREREE ERE KRREEERKKERERREKKRERKERKEKKAAKK ee 


* 


* 


SUBROUTINE LOAD85 


SUBROUTINE TO DEVELOP A LOADING, OR VORTEX DISTRIBUTION, THAT 
IS ASSOCIATED WITH A FLAT PLATE WING. 
IS RELATED TO THE PRESSURE DIFFERENTIAL, DELTA CP, AND HAS THE 
SAME BASIC SHAPE. 


OF INFLUENCE COEFFICIENTS AND THE WING SHAPE, AND SOLVING FOR 


THE VORTEX STRENGTHS. 


ALL REQUIRED PARAMETERS AND VARIABLES ARE CONTAINED IN COMMON. 


A DESCRIPTION OF THE COMMON BLOCK VARIABLES IS CONTAINED IN 


OTHER ROUTINES. 


THE SOLUTION IS GAINED BY USING A MATRIX 


VARIABLES USED IN THE SUBROUTINE ARE: 


ALF 


DELTA 


ADELTA 


CBAR 


CL 


CD 


CM 


CXJ 


CZJ 


ANGLE OF ATTACK IN RADIANS 


THE LEADING EDGE SWEEP 
THE LEADING EDGE SWEEP 
THE AVERAGE CHORD 

INTERIM VALUE, USED TO 
INTERIM VALUE, USED TO 
INTERIM VALUE, USED TO 
INTERIM VALUE, USED TO 


INTERIM VALUE, USED TO 
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ANGLE IN RADIANS 


ANGLE IN DEGREES 


COMPUTE A FINAL VALUE 


COMPUTE A FINAL VALUE 


COMPUTE A FINAL VALUE 


COMPUTE A FINAL VALUE 


COMPUTE A FINAL VALUE 


THE LOADING GENERATED 


* XAC 

* CLT 

* ayn 

* CMT 

* 

* CDOCL2 
x 

* SUBROUTINES 
* SET85 
* 

* DNWS85 
x* 

* 

* GAUSS 
* 

x 

- MULT2 


INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 


: X POSITION OF THE AERODYNAMIC CENTER 


PeletAte Gur COEFFICIENT OF THE FLAT WING 


TOTAL DRAG COEFFICIENT OF THE FLAT WING 
TOTAL MOMENT COEFFICIENT OF THE FLAT WING 
ABOUT THE AERODYNAMIC CENTER 


GO7GL°2, ALSO CALLED K IN THE DRAG POLAR 


CALLED ARE: 

FORTRAN, USED TO GENERATE THE WING AND SECTION 
ARRAY CONTAINING PLANFORM COORDINATES. 

FORTRAN, USED TO COMPUTE THE DOWNWASH GENERATED 
AT A POINT BY A HORSE-SHOE VORTEX AND IT'S MIRROR 
PMAGESON THE OPPOSITE SEMI-SPAN. 

FORTRAN, USED TO SOLVE THE INFLUENCE COEFF MATRIX 
AND SLOPE VECTOR, FOR THE VORTEX STRENGTH. USES 
SIMPLE GAUSS ELIMINATION. 

FORTRAN, USED TO COMPUTE THE DOWNWASH GENERATED 
ON EACH VORTEX ELEMENT, NECESSARY TO FIND THE 


FORCES ACTING ON THE WING. 


* 


FHKE KAKKKKK KEK KEK ERK KEK KEKE KEE BAKKE EEE EK EEG BK KEK HHH HHH HHH EHH HHH KE KE EH 


SUBROUTINE LOAD85 


REAL LAMDA 


CHARACTER*3 , 


ELLIP 


COMMON PI,B,WING(101,23), SECTN(11,10) ,AR, LAMDA, DELTA, NX,NY, 


fae 


* ALPHA,CLDSRD,ELLIP,XAC,CLA 
COMMON /COF/A(101,101), NEQNS 
COMMON /CAM/AA( 101,101), AS(101,101) 
COMMON/FILS/IIN, IOUT ,JOUT , KOUT , LOUT , MOUT , NOUT 


60 FORMAT ( ’ FLAT WING’, 


A/T = ' F12.6/'’ CD = ',F12.6/' CD/CL2 = ',F12.6 
* /' CMAC = ',F12.6/' XAC = ',F12.6/’ AR = ', F126 
* /' LAMDA = ',F12,6/’ DELTA = "(Flan ae’ 13 

= /' NY att ont /' ALPHA = ‘,F12.06/' €LDSRD — 3 r lee 
x  /' ELLIP = ',A3 #£x/' CLa, Lift Curve Slope’ 


* /' per/Deg= ‘',F12.6) 


KHKKK KKK KKK KEK KK KK KKK KK AK KEK RE KKK KERK KARE KKK KEKE KEK AK KEK KEK BERK E KER EEE 


* * 
* COMPUTE THE ANGLE OF ATTACK IN RADIANS. * 
* * 


KREKKKRKREKEKK KKK KE KEKE KK KE KEKKKKHK HK HK HK HHH KH K Hk Kk kk ek HK HK ok vk ok ok Kk ok kk Kk Kk kK 
ALF = ALPHA*PI/180.0 


SINALF 


SIN(AEF) 


COSALF 


COS(ALF) 


KKK KKK KKK KKK KKK KKK KKK KEK KK KK ERE KK KKK KEKE KKK EKER RE KEKE KEK KE KE ER RE KEK KEK KKK 


* * 
* CALL THE SETUP ROUTINE. USED TO BUILD THE COORDINATES OF THE * 
* WING. THEY ARE STORED IN THE ARRAYS WING AND SECTN. * 
* * 


IGG 


KKK KK HHH HHH HH WKH Hea HEH ae EK aE EE aE ak Ok BE IK ak HK ak vk ME ak ak Oe OE ae ak ah ak ak ak Oh BE De EE AE AE IEE EME IE EE ME BEDE DE KE aE MEE KE 


CALL SET85 


Wak akeok kako kok ak ak ak ak ok ek ak ok ae & sh ak ok ak ok a ak ok ak ok ak ak Oe ak ek ak AE ak ak EE HEM ae kM ak EK ae EE BEE ek EH a EH I EEE 


* * 
* BUILD THE THREE COEFFICIENT MATRICES, A(), AA(), AND AS() * 
* A() AND AA() ARE THE SAME, AND ARE FOR EVALUATING FLOW * 
* TANGENCY AT THE EDGE OF THE GRID ELEMENTS. * 
* AS() IS FOR EVALUATING DOWN-WASH GENERATED ON THE VORTEX * 
* ELEMENT ITSELF. 

* * 


KREKAK KKK KEE HKK KKK KKH KHL K KH Kw Hk ak sek se sk ak kk ak vk ok sek ak Hk sek vk sk Hh ak ak ok sk vk shee se he se eK K ake ke ek ke 
DO 130 I=1,NY 
DO 120 J=1,NX 
TJ = (I-1)4*NX4+J 


HH WKHHR AKA AEH NA KK KAN AK KKH HHH KKK HH HHH HH HH HHH HH HH HH HHH a-ak vk 9 ae ok se ak ak ak vk ok ok ak ak vk oh ok 


* * 
* AUGMENT THE A() MATRIX WITH THE RIGHT HAND SIDE, WING SHAPE. * 
* * 


HRAAHHAH HANH AH WH WH HH HHH HHH HHH HHH H WH WH WW Wk HW HK sek ak ae ak HK ak ak ok ak vk ok ak ok ak vk ok vk ak ok ok ok ok ak ok ok ok ok 


A(IJ , NEQNS+1) 


I 


SINALF 


* A(CIJ ,NEQNS+1) ALF 


DO 110 K=1,NY 
DO 100 L=1,NX 
KL = (K-1)*NX+L 


CALL DNWS85(1J,KL,AC(KL,IJ),*CONT’ ) 


abey) 


AACKL,IJ) = A(KL,IJ) 


CALL DNWS85(IJ,KL,AS(KL,IJ),‘SELF’) 


100 CONTINUE 
110 CONTINUE 
720 CONTINUE 


VEEN, CONTINUE 


KKHK KKK KK ak Hk KKK KKK HK HK HHA HAHA HHA AH HH Hk ok ok Hk ok ak ok ak ak vk HM Mk ak vk vk vk Kk HH HK HK AK kk kk vk kk KK 


* * 
* COMPUTE THE AVERAGE CHORD * 
* * 


KHREKEKAKEK KEK KEK KK EKER KKK KKK KK KKK KKK EK KKK KEK KKK KK EK KKK KKK KK EK KKK DEK KK EK ERY 


CBAR = (SECTN(1,1)+SECTN(NY,1))/2.0 


KHAKKKKKK KKK KAKA A AAA KKK KAKA AK AKAK KK KA KKK KK AKK KKK KAKA KK KK KK KE KKK KK KAK KK KAKE KEEN 


* * 
* CALL THE GAUSS ROUTINE, TO SOLVE FOR THE LOAD VECTOR * 
* THE RESULT IS RETURNED IN PLACE OF THE ORIGINAL LOAD VECTOR * 
* * 


HKHKAHKHH KAKA KHAN AK KKK KK AHA AKA KK HK HH KHH HHH HH HH HHH HHH HHH HH HH HHH HHH HH KWH K KK F 


CALL GAUSS(1) 


KARE KAEK KKK AKER KEK KKK KKK KKK HK KKK KKK HH KH HHH HH HK HW HK KWH HV KKK AH HAAN KALA AK KAA AAS 


* * 
* STORE THE RESULT IN THE WING() MATRIX * 
* CONVERT THE CIRCULATION, OR VORTEX STRENGTH IN ELEMENT 15 TO * 
* DELTA CP, IN ELEMENT 16 x 
* * 


KAHHK HH HHH HAH HHH HHH HW HHH HHH HH WHHHH HHH HHH HH HH WH HHH HHH RH Hk HHH K Hk HHH HK HK KK KA K 


‘oe 


&=> 


Dor 10 (>= l,.NY 
DO 20079 1, NX 


TJ = (1-1) 4Nx+J 


WING(IJ,15) = A(IJ,NEQNS+1) 


WING(1IJ,16) 


2.0*WING(IJ,15)/WING(IJ , 13) 
200 CONTINUE 
210 CONTINUE 


KRKRAKKKH KEK KKH A HAHAH AKA KHANH HK AK AKA K HHH HK HHH HK HAHAHA KAHN AKAN A AK KAA KAKA NK AAA EK 


* * 
* COMPUTE THE DOWNWASH ON EACH VORTEX ELEMENT, NECESSARY TO * 
* PINDeTHE FORCES ACTING ON THE WING. * 
* * 


KHEAAAAA KAKA AKA AKAN KAA KA KAKA AK KKH KKK KKH AHHH AHHH NAH HK KK KAKA AA HK KKK KKK HAA AA KEK 


GALL, MULT2 


KAEKAAAAAKEAA AEA A AKA EAA AERA AER ERAN KEKE KH KEK KEHNA ANAK AK GK NHN AAG A KK AAA A EEK 


* * 
* INITIALIZE VALUES AND COMPUTE THE FORCES ACTING ON THE WING * 
* SECTIONS AND ALSO THE TOTAL WING. * 
* * 


HARKER KKAK EAA AK EAA EAN KEAKAKK EEK EK KEK KEE KEKE KKH KHER HHH KKH HHH HHH HE HH HK HK MH ee 


CD = 0.0 
CL = 0.0 
CM = 0.0 


DO 330 I=1,NY 
CX¥J = 0.0 


CZJ = 0.0 


el 9 


CMJ = 0.0 
DO 320 J=1,NX 
IJ = (1-1) *NX4J 


CZ 


WING( IJ ,15)*2.0*(1.-WING(IJ ,22)*SINALF) /(B*CBAR) 


CXJ 


I 


WING(IJ,15)*WING(IJ , 22) *2 .O*COSALF/ (B*CBAR) 


CMJ 


I 


-WING(IJ,15)*WING(IJ ,4)*2.0* 


* (1.-WING(IJ ,22)*SINALF) / (B*CBAR*CBAR) 


SECTN(I, 2) SECTIN(1,2)+CZJ 


SECTN(I,3) SECTN(I,3)+CXJ 


SECTN(I,4) 


SECIN(I ,4)+CMJ 


520 CONTINUE 


CL = CL+SECTN(I,2)*SECTN(I,6) 


CD CD+SECTN(1,3)*SECIN(1,6) 


CM 


CM+SECTIN(1I,4)*SECTN(1I, 6) 
SECIN(I,5) = -SECTN(1,4)/SECIN(1,2) 
330 CONTINUE 


HHKKEKAK HHH KKH KH HHH HH HH HMA AA AK AK HHA KA HAH KKH HH HHH KAKA KKH HA HH HK KK KK KH KA KKK HK 


* * 
* NOW COMPUTE THE MOMENTS ABOUT THE AERODYNAMIC CENTER * 
* * 


KKEKKKHK HHH HHH HHH HH HHH KKH EHH K AK KKH K HEE H KK KH HEHEHE KK KKE KK KE KK KAKA KH KKK KE KEK 
XAC = -CM*CBAR/CL 
DO-335 = Ny 
SECTN(1,4) = SECTN(1,4) + SECTN(I,2)*XAC/CBAR 
335 CONTINUE 


KK KK KKH HHH HHH HHH KKH KKK HH HK KH HH HHH HH HHH HH HH HH HHH HK HH HH HH HH HHH HH AHHH HH HH WK 


120 


* COMPUTE THE FINAL TOTAL WING COEFFICIENTS AND WRITE TO THE * 
* 1Q DE} S * 
* * 


KKEKKEKKE KARE KEKE KEK AK KAKEKK EK KEKE KEKE KE RE RAK KEKE KARA ERE RK EER ERK EKER AE KEK 


Cri wOxCL 


CDT = 2.0*CD 


CMT 


2.0*CM + (CLT*XAC/CBAR) 


CLA = CLT/ALPHA 


CDOCL2 


CDT /CLT**2 


ADELTA 


DELTA * 180./PI 
OPEN (UNIT = MOUT, FILE ='FLAT.DAT’) 
WRITE(MOUT,60) CLT,CDT,CDOCL2,CMT,XAC,AR, LAMDA, ADELTA, NX, NY, 
* ALPHA, CLDSRD, ELLIP, CLA 
CLOSE (MOUT) 
340 RETURN 


END 


ileal 


SHAP85 


KKKKKKKKKKK KKK KEKE KK KEKE K KKK KEKE KKKEKEKKKEKKKEK EKER ERE KK KKK KK KEKE KE KEKEKKEK 


* 


* 


SUBROUTINE SHAP85 


SUBROUTINE TO DEVELOP A WING SHAPE THAT IS ASSOCIATED WITH 

A SPECIFIED LOADING FUNCTION. THE LOADING FUNCTION IS USUALLY 
ELLIPTIC, IN THE SPAN AND CHORD DIRECTION. THE SHAPE Wite 
USUALLY RESEMBLE A PARABOLIC CAMBER, IN CHORD, AND HAVE SOME 


DEGREE OF TWIST, ALONG THE SPAN 


ALL REQUIRED PARAMETERS AND VARIABLES ARE CONTAINED IN COMMON. 
A DESCRIPTION OF THE COMMON BLOCK VARIABLES IS CONTAINED IN 


OTHER ROUTINES. 


VARTABLES USED IN THE SUBRGUTINE TARE: 

CBAR : THE AVERAGE CHORD 

CLILST “> GLREF FROM THE PREVIOUS ZITERAIEG 

CLREr : A SCALER, USED TO INCREASE OR DECREASE THE AMPL 
OF THE SHAPE, IN ORDER TO ACHIEVE THE DESIRED 
LIF? -COEFLICUE ae 

MOD : A COUNTER USED TO LIMIT THE NUMBER OF ITERATIONS 
THAT THE PROGRAM CAN USE IN AN ATTEMPT TO GENERATE 
A DESIRED LIFT COEFEICrENa 

H : A TEMPORARY VARIABLE USED WHILE INTEGRATING THE 


HEIGHT OF THE WING ALONG A CHORD. 


ZZ 


* 


* 


* 


* 


% 


% 


% 


* COSCAM : 


* SINCAM : 
* 

* CCL 

* CCD 

* CCM 

* CXCP 

* CCL. 

* CCDT 

* CCMT 

* 

* CCDOCL : 
* 

* SUBROUTINES 
* MULT1 

* 

* 

* MULT3 

* 

* 


THE COSINE OF THE ANGLE MADE BY AN INDIVIDUAL GRID 
ELEMENT AND THE REMOTE VELOCITY. 

THE SIN OF THE ANGLE MADE BY AN INDIVIDUAL GRID 
ELEMENT AND THE REMOTE VELOCITY. 

INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 
INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 


INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 


;: X POSITION OF THE CENTER OF PRESSURE 


TOTAL LIFT COEFFICIENT OF THE CAMBERED WING. 
TOTAL DRAG COEFFICIENT OF THE CAMBERED WING. 
TOTAL MOMENT COEFFICIENT OF THE CAMBERED WING. 
ABOUT THE AERODYNAMIC CENTER 


CD/CL*2, ALSO CALLED K IN THE DRAG POLAR 


CALLED ARE: 

FORTRAN, USED TO COMPUTE THE SHAPE OF THE WING 
FROM THE LOADING VECTOR AND THE INFLUENCE COEFF 
MATRIX 

FORTRAN, USED TO COMPUTE THE DOWNWASH GENERATED 
ON EACH VORTEX ELEMENT, NECESSARY TO FIND THE 


FORCES ACTING ON THE WING. 


* 


* 


HHH WH KAHA AHH HHH HHH HHA AHA WHA A KAA WAN HK K HANH AHH AHA AHA AA KKK KEKE KAKA KKK KEK KEENE 


SUBROUTINE SHAP85 


REAL LAMDA 


CHARACTER*3, 


ELLIP 


eZ 3 


66 


70 


COMMON PI,B,WING(101,23), SECTN(11,10) ,AR, LAMDA, DELTA,NX,NY, 


* ALPHA, CLDSRD,ELLIP,XAC,CLA 


* 


5 


COMMON/COF/A(101,101), NEQNS 


COMMON/CAM/AA(101,101), AS(101,101) 


COMMON /FILS/IIN, LOUT,JOUT ,KOUT CLOUT Metter 


CEarESt 


CLREF = 0.05 


MOD = 


= 0.0 


0 


FORMAT (’ UNABLE TO REACH THAT CLDSRD’ ) 


FORMAT ( ’ 


: 
, 
2 
* 
- 


CL 


CMAC 


NY 


SA ih 


, 


$ 


, 


, 


F12.6/! 


1267, 


Bats 


aS 


,A3) 


Vis 


CAMBERED WING’, 


cD = ' P1267’ CD/CL2 =F lome 
' AR = ' Plgee 

DELTA = ',F12.6/' NX = ' 13 

ALPHA = ',F12.6/' CLDSRD = ',F12.6 


HHH AH KHK KKK HAH WK WH HH HH HK Vv Wk a vk vk ok a vk HK ak vk ok 9k vk vk ME 9k 9K ak Dk a ak vk Vk VE IE IK ak DE Ok OH ok ak OE OK ok KK 


* 


* 


* 


COMPUTE THE AVERAGE CHORD 


* 


* 


* 


PKK KK KE EK HK IE HE EK aK ak vk EH ak ak I EAE OK aE aE IK OE 2k ak De Ob ak ok ak ob ob 2k ak ob ak ak ak ak ob ob ok 2k ok oh ok ak ak ak ok ok ak ak oh oh ok ak ok Oe OE a Ok EK 


CBAR = (SECTIN(1,1)+SECTN(NY,1))/2.0 


EK HH HHH HK KE HHH HEE HE EE EME AE IE aE IK ah ak I 2k oe ak oh ah 9k 3 ak 3h ah 9k 3h 9k 3] 9h ak 3b ak ok ob 2k 9 ak ok ak ok ah ak oF ak of ak ok ak of ok ok OF HK OK 


* 


* 


USE THE SCALER, CLREF, TO PUT AN ELLIPTIC LOAD DISTRIBUTION 


IN THE VECTOR (IJ,17). 


IN THE VECTOR (IJ,18). 


AT THE SAME TIME COMPUTE THE DELTA CP 


line 


* 


* 


a 


* 


* 


KKK KKK KKK KKK KH KKK KEKE KEK KEKE KEKE KEK KKK AK KEKE KEK EK KK EK EKER AER KKK RRR ERK ERE RK 


390 


400 


DO 400 IJ = 1,NEQNS 


WING(IJ,1/7) WING(IJ ,21)*CLREF 


WING(IJ,18) 


2 .O*WING(IJ,17)/WING(IJ,13) 


CONTINUE 


KKKKIK KK KEK KKK KEKE KEKE KEKE KEKE KEK RE KEK KER ER KE RK ERE ER EER KEKE RR ERE EERE KEKE 


* 


* 


* 


* 


* 


CALL THE SUBROUTINE MULT1 WHICH WILL 
MUL Tiebweehe LOAD VECTOR, (1J,17) BY THE COEFFICIENT MATRIX, 


AAC ye) TO GET THE SHAPE VECTOR, (1J,19) 


* 


* 


* 


* 


* 


HEKKEKKKEKEKEKEKEK KEKE KEKE KEKE KKKEKEK EKER KERKE ERE KKK REE RE RKKE RK KEK ERE KERR REESE 


CALL MULT 


KEEKKKKKKEKAE EE ERK ERE KKERK EKER ERK KERR ER ER ERE RRE RRR RERRREREREERERRERERERREREEESD 


* 


* 


* 


* 


CHECK THE RESULTING SHAPE VECTOR (1J,19), TO SEE IF THE ANGLE 
THAT WAS REQUIRED HAS A SIN GREATER THAN 1, THAT IS SIMPLY 
CHECKING THE SHAPE VECTOR, TO SEE IF IT IS GREATER THAN 1. 

IF IT IS, THE SCALER MULTIPLIER, CLREF, IS DECREASED. AT THE 
SAME TIME, MOD, THE ITERATION LIMITER, IS INCREMENTED. IF 
MOD IS GREATER THAN 5, THAT MEANS WE OVERSHOT 5 TIMES, AND 
PROBABLY WILL NEVER GET TO THE DESIRED LIFT COEFFICIENT. 

IN THAT CASE, COMPUTE THE LAST VALID SHAPE, AND PROVIDE THOSE 


DATA. 


125 


* 


Pag 


e 


KKK RHK RRR HK KHHH KKH KK HHKHHHW HHH HH RRR HR WKH K HH HHH RRR KKH HHWWHKHKA KKKKKKKA KKKKKK 
DO 410 IJ = 1,NEQNS 
IF (ABS(WING(IJ,19)).GT.1.0) THEN 
CLREF = CLREF*.9 
MOD = MOD + 1 
IF (MOD.GT.5) THEN 
WRITE (IOUT, 66) 
CALL DELAY (6) 
GO-TO s5u0 
ELSE 
GO TO 390 
ENDIF 
ELSE 
ENDIF 
410 CONTINUE 
HHKKKE KAKA KEKE KEKE KKK KEK KERR KEKE RK KEKE KK HE KKK KKH KH KHER KK HEE KEKE KEKE KS 


* 


% 


& BEGIN TO COMPUTE THE FORCE AND MOMENT COEFFICIENTS ON THE * 
‘ CAMBERED WING. MULT3 WILL MULTIPLY THE SPECIFIED LOADING * 
* BY THE COEFFICIENT MATRIX, AS, TO GET THE DOWNWASH AT THE * 
* CENTER OF EACH CAMBERED ELEMENT. INTERMEDIATE VARIABLES ARE * 
* INITIALIZED. : 
* COEFFICIENTS FOR SECTIONS AND THE WHOLE WING ARE COMPUTED. * 


* 


% 


KKAKAN AKAN AK AK KKK AKER AK KK KB RKE KEK KKK KKK KE KEKE KEKE KE KK ER KK NK KE KER KEKE RK ENE EK 


CALL MULT3 


126 


20 


CoDe= O70 
CCL = 0.0 
CCM = 0.0 


DO 530 I=1,NY 


CCXJ = 0.0 

CCZJ = 0.0 

CCMJ = 0.0 
SECTN(I,7) = 0.0 
SECTN(1,8) = 0.0 
SECTN(1,9) = 0.0 
H = 0.0 


IJ = (1-1)*NX+1 
DO 520 J=1,NX 
IJ = (1-1) *NX+J 
H = H-WING(IJ,19)*WING(IJ,13) 
WING(IJ,20) =H 
COSCAM = COS(ASIN(WING(IJ,19))) 


SINCAM 


WING(IJ ,19) 
CCZJ = WING(IJ,1/7)*2.0*(1.-WING(IJ ,23)*SINCAM) /(B*CBAR) 
CCXJ = WING(IJ,1/7)*WING(IJ ,23)*2.0*COSCAM/(B*CBAR) 


CCMJ = -WING(IJ,17)*WING(IJ ,4)*2.0* 


* (1.-WING(IJ ,23)*SINCAM) /(B*CBAR*CBAR) 


SECUN( tay jem SECIN( 157 )+CCZJ 


SECTN(1,8) = SECTN(I,8)+CCXJ 
SECTN(I,9) = SECTN(1I,9)+CCMJ 
CONTINUE 


ay 


CCL = CCGL+SECIN(1 (7) *SECINGI 
CCD = CCD+SECTN(1I,8)*SECTN(I, 6) 
CCM = CCM+SECTN(I,9)*SECTN(I, 6) 


KEKE KAKKAK KK KAK KEK KE KKK KKK KE KEKE RE EEK EKER KEE KK KER ERR ERE RERK KEKE RE RE REE RER 


* * 
* COMPUTE THE SECTION TWIST. THE CAMBERED WING GE * 
* WILL BE THE SAME AS CLDSRD, WITHIN A VERY SMALL MARGIN. * 
* * 


HAEKKAAKEKEK KEK KKK KEK KEK KEKE EKER EKER KEK KEKE KE KEKE KER KEKE RAE KRER ERE REE KERR 


SECTN(I,10) = -WINGCI*NX, 20)*180./(SECTN(I,1)*PI) 


HHA AEKEKRERE RE RA KAKA AKER KEKE KE KEKE KKK KEK KKEKEKKKKKERKERKERREERREEKEKEKEEEKERER 


* * 
* COMPUTE THE MOMENT COEFFICIENT ABOUT THE AERODYNAMIC CENTER. * 
* 7 


HKKA KKK HAHAHA KKK A KAKA AA KK KKH K KKK KKK ER DER HHH TER ER KW WKH HHH HK 
SECTN(I,9) = SECIN(I,9) + (SECTN(I, 7)*XAC/CBAR) 
530 CONTINUE 


CXCP = -CCM/CCL 


CCLT 2.0*CCL 


CCDT 2.0*CCD 


CCMT 


2.Q0*CCM + (CCLT*XAC/CBAR) 


CCDOCL 


GCDI7 Gel T**2 


GLilsT CLREF 


IF (ABS(CLDSRD-CCLT).LT.(CLDSRD*0.01)) GO TO 540 
CLREF = CLREF*CLDSRD/CCLT 


GO TO 390 


128 


FEHR AHA AHHH HAHAH HK KH KHAHHAHA HHA KH K HAH KH HH HK KH KK KKK KKK KKK KK KEK EK KEK KEK EERE 


* * 
* THE DESIRED LIFT COEFFICIENT HAS BEEN ACHIEVED, WITHIN A * 
* A TOLERANCE OF +/- 1.0%. THE CORRECT VORTEX STRENGTH AND * 
* DELTA CP ARE PUT INTO THE WING ARRAY, AND THE FINAL CAMBERED * 
* WING VALUED OF LIFT, DRAG, K, MOMENT, AND CENTER OF PRESSURE * 
* ARE PRINTED. * 
* * 


FRI KKA KAKA KAA AKA KAKA KE KAKA RK KKH HHL HHH KH HEE EHH HEE KE HK HEE 
540 CONTINUE 


DO 550 IJ = 1,NEQNS 


WING(1IJ,17) = WING(IJ,21)*CLILST 


WING(IJ,18) = 2.0*WING(IJ,17)/WING(IJ,13) 
550 CONTINUE 


KRHKKANKEKEKKAKEK AK KKK KEK KKK EK KKK KEK EKER EEK EKER EKER ER ERK EEE REE EK KEKE EEE EEE 


* * 
* WRITE THE TOTAL WING COEFFICIENTS TO THE SCREEN * 
* * 


KKEKRAKKAKAKKEKEKK KKK AK KEKE KEK KARE RE i bb bbb bbb ik bsb ot bk ik bb bok bot ok bab ok bak ok 
OPEN (UNIT = MOUT, FILE = 'CAMBER.DAT’) 
ADELTA = DELTA *180./PI 
WRITE(MOUT, 70) CCLT,CCDT,CCDOCL, CCMT, AR, LAMDA, ADELTA, NX, NY, 
* ALPHA, CLDSRD, ELLIP 
CLOSE (MOUT) 
RETURN 


END 


129 


SET85 


KAEKREKREKKEKE KEKE EKER ERR EKER KKK RRR A RRA KERR EKER RRR REAR KAEKKRKRE KEES 


* * 
* SUBROUTINE SET85 * 
* * 
* This program is written for the setup of a straighthorse shoe * 
* vortex over tapered and swept wing. It uses cosine spacing * 
* for both the span and chord elements. This concentrates * 
* elements near the edges of the wing, where they are most * 
* important. * 
* * 
* Column elements in the WING array are: * 
* l sequential position, also called IJ * 
* 2 integer y position 2 
* 3 integer x position * 
* 4 x center of element * 
* 5 y center of element * 
* 6 % position of horse shoe vortex * 
* 7 blank * 
* 8 ya position of one corner of horse shoe vortex * 
* 9 yb position of one corner of horse shoe vortex * 
* 10 xp position of flow tangency point, 3/4 chord * 
* ll xp position of down wash point, 1/4 chord * 
* 12 yp position of flow tangency or down wash pt, center of elem * 
* 13 del x, chord wise dimension * 
* 14 del y, span wise dimension > 


hs 


16 


a 


te 


In, 


20 


2): 


ZZ 


Z3 


Elae plate vorticity 

flat plate delta cp 

cambered voricity, adjusted for CLi 
cambered delta cp 

cambered slope 

cambered height 

cambered vorticity, if CLREF were = 1.0 
downwash at center of flat plate element 


downwash at center of cambered element 


Column elements in the SECTN array are: 


1 


iz 


10 


chord 

flat wing lift force per unit span 

flat wing drag force per unit span 

flat wing moment about the aerodynamic center 
flat wing XCP, center of pressure relative to x=0 
section delta y 

cambered wing lift force per unit span 

cambered wing drag force per unit span 

cambered wing moment about the aerodynamic center 


cambered wing twist. 


VARIABLES USED ARE: 


B >: WING SPAN 
DTHE > AN ANGULAR ELEMENT WIDTH, IN RADIANS ' SPAN 
DPSI > AN ANGULAR ELEMENT DEPTH, IN RADIANS, CHORD 


ema 


x 


% 


* THEA > INTERMEDIATE VARIABLE * 


* THEA : INTERMEDIATE VARIABLE * 
* PSIA : INTERMEDIATE VARIABLE * 
* PSIB : INTERMEDIATE VARIABLE * 
* ELE : INTERMEDIATE VARIABLE * 
* ELT : INTERMEDIATE VARIABLE * 
* PSI : CHORDWISE ANGULAR COORD FOR THE CENTER OF AN * 
* ELEMENT. USED IN GENERATING ELLIPTIC LOAD. * 
* THETA : SPANWISE ANGULAR COORD FOR THE CENTER OF AN * 
* ELEMENT. USED IN GENERATING ELLIPTIC LOAD. * 
rs x 
* FUNCTIONS USED ARE: 

* BETA : COMPUTES THE LENGTH OF A LOCAL SECTION CHORD. * 
* * 
* ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON BLOCKS * 
* THE ONLY VARIABLES CHANGED IN THIS ROUTINE ARE: * 
*« B * 
* WING() * 
* SECTN() * 
re * 


KA KH AKA AK KAA AHH AK KKK A KK ANAK KA AK HE AKA A HEA HELA HEA EAA AEA KEK AK HHA HR HN HK 
SUBROUTINE SET85 
COMMON PI,B,WING(101,23), SECTN(11,10) ,AR,LAMDA, DELTA,NX,NY, 
* ALPHA,CLDSRD,ELLIP,XAC,CLA 
REAL LAMDA 


KRKKA KEE KK KKK AK EK KEELE KK KKK EEK EEK KKK ELEN EAA K EEK EEK EEK KEKE KK EAN SE AEE AY 


OZ 


* DEFINE THE FUNCTION BETA, WHICH COMPUTES THE LENGTH OF A LOCAL * 
* CHORD. * 
* * 


KKK KA KEKE EEK KEKE KEK KEE KEKE ERE ERRRE KERR REE RERERERERRRERR REE ENR 


BETA(Y,LAMDA,AR) = 1.0-(Y*4.0*(1.0-LAMDA) /(AR*(1.0+LAMDA) ) ) 


KAKKKKEEKKKEKKEKKEKEKEKEEEKEREER KEKE RREKREERKEEERKRKEEK KERR KRERE KEKE REE KEREEREREERERE 


* * 
* DEFINE THE SIZE OF SPAN AND CHORD WISE GRID ELEMENTS. * 
* DETERMINE THE TOTAL WING SPAN. * 
* * 


HHH AKA KKKAAK AKA AERA KAKA ERER 
* DTHE = PI*0.5/FLOAT (NY) 


KEKEEKEKEEKEKEKEE KEKE KERR KEKRERERERERERKRERKRKRERKKREKREEREEREKRERERRRRRRERKKRK 


x * 
* THIS STATEMENT SETS THE GRID UP SO THAT THE OUTER ELEMENT * 
* IS IGNORED. i 
* * 


KKEKAKEKAEKEKREKERERE EERE KEKE RRER RRR ERR RRR RRR RERRERER RRR ERERERER 


DTHE PI*O.5/FLOAT (NY+1) 


DPSI = PI/FLOAT(NX) 
B = AR*(1.0+LAMDA) /2.0 
DO 110 I = 1,NY 

DO 100 J = 1,NX 


KKKEKKERKEKKEKEAKK KK KK KEKE KEKE EKER KEE KRER EER REE EKER ER RERERREKRERRERER EERE 


* 7 


33 


* COMPUTE ALL REQUIRED COORDINATES FOR THE WING() ARRAY. * 
* * 
KKKA AAA AAA HAHAHA HHH HAHA KAA HAHAHAHAHA KKK KKK A KAHN HK HK KH KH HA KHAN KKK 


IJ = (I-1)*NX+J 


WING(IJ,1) = FLOAT(IJ) 

WING(IJ,2) = FLOAT(I) 

WING(IJ,3) = FLOAT(J) 

THEA = (FLOAT(I-1)*DTHE)+(PI/2.0) 


THEB 


(FLOAT(1I)*DTHE)+(PI/2.0) 


PSIA 


(FLOAT (J-1)*DPST) 


PSIB (FLOAT(J)*DPSTI) 


WING (IJ ,8) -B*COS (THEA) /2 .0 


WING(IJ,9) 


-B*COS (THEB) /2.0 
WING(IJ,12) = -B*(COS(THEA)+COS(THEB) ) /4.0 
ELE = (1.0 - COS(PSIA))* 
BETA(WING(IJ,12),LAMDA,AR)/2.0 
Bish =o 0 = COS CES 1B) > 
* BETA(WING(IJ,12) , LAMDA, AR) /2.0 
* THESE TWO STATEMENTS SET THE POINTS AT 1/4 AND 3/4 CHORD 


WING (IJ ,10) 


WING (IJ ,12)*TAN (DELTA) +(ELE+3*ELT) /4.0 


WING(IJ,11) 


WING(IJ ,12)*TAN(DELTA)+(3*ELE+ELT) /4.0 


WING(IJ,4) = WING(IJ,12)*TAN(DELTA)+(ELE+ELT) /2 .0 
WING(IJ,6) = WING(IJ,11) 
WING(IJ,5) = WING(IJ,12) 


WING(IJ,13) ELT-ELE 


WING(IJ,14) = WING(Is, 9) -WiniaGmamer 


134 


PSI = ACOS(1.0-(2.0*(WING(IJ ,11)-(WING(IJ ,5)* 
* TAN(DELTA) )) /BETA(WING(IJ ,5) , LAMDA, AR) ) ) 

THETA = ACOS(-2 .O*WING(IJ,5)/B) 

WING(IJ,21) = 4.0*(1.0+LAMDA) *SIN( THETA) *SIN(PSI)* 
* WING(IJ ,13)/(PI**2*BETA(WING(IJ,5) , LAMDA, AR) ) 


WING(IJ,7) = 0.0 


| 
© 
=, 


WING(IJ ,15) 
WING(IJ,16) = 0.0 
WING(IJ,17) = 0.0 
WING(IJ,18) = 0.0 
WING(IJ,19) = 0.0 
WING(IJ,20) = 0.0 
WING(IJ,22) = 0.0 
Wane Ci J. 25) = 020 
100 CONTINUE 


RREKKREKER AKER EKA ER ARR RE REE RR EAR ERE BSE EKA EERE RAN BEEK E RENE KKK K KEE AK KEN 


% 
nN 
x 
- 


* * 
* COMPUTE ALL REQUIRED COORDINATES FOR THE SECTN() ARRAY. * 
* * 


BRC CI SK IK KARR KR RRR RK RK RR RRR ER EREKRERRERRRERERE 


SECTN(I,1) = BETA(WING(IJ,12),LAMDA,AR) 
SECTIN(I,2) = 0.0 
SECTN(I,3) = 0.0 
SECTN(I,4) = 0.0 
SECTN(I,5) = 0.0 


SECIN(1,6) = WING(IJ,14) 


3:5 


SECTN(I,7) = 0.0 
SECTN(I,8) = 0.0 
SECTN(I,9) = 0.0 


SECTN(I,10) = 0.0 
110 CONTINUE 
RETURN 


END 


Ls 


DNWS85 


FREAK KEKE EE KK KK KKK KK ERK ERE RE KE KEK ERE KE KKK EK RR ER RE KK EKER ERR ERE RR ER ER 


* 


* 


SUBROUTINE DNWS85(IJ,KL,W, IND) 


COMPUTE DOWNWASH ON GRID ELEMENT NUMBER KL DUE TO 


VORTICES ON GRID ELEMENT IJ, AND IT’'S MIRROR IMAGE 


VARIABLES USED ARE: 


IJ ; VORTEX GRID ELEMENT NUMBER, INPUT 

KL > DOWNWASH GRID ELEMENT NUMBER, INPUT 

W : THE DOWN WASH, RETURNED TO THE CALLING PROGRAM 
OUTPUT 

IND > A FLAG TO DETERMINE WHETHER TO USE TRAILING 


EDGE OF THE GRID ELEMENT, OR THE CENTER OF 


DHE SECEMENT (> INPUT 


FUNCTIONS USED ARE: 
WHV1 : FORTRAN, COMPUTES THE DOWNWASH DUE TO ONE CORNER 
OF THE HORSESHOE VORTEX 
WHV2 : FORTRAN, COMPUTES THE DOWNWASH DUE TO JUST THE 


TRAILING FILAMENTS OF THE HORSESHOE VORTEX 


* 


* 


* 


* 


* 


WH KKK KKK KEKE KEKE KEE RK ERE RRA K KEK EERE KR ER ERE RE R RRR RRR RE REE RRR EEE EEE EE EEE 


SUBROUTINE DNWS85(IJ,KL,W, IND) 


Jay 


COMMON PI,B,WING(101,23), SECTN(11,10),AR,LAMDA, DELTA, NX,NY, 
* ALPHA, CLDSRD,ELLIP,XAC,CLA 

CHARACTER*4 IND 

IF(IND.EQ.’CONT’) GOTO 100 


EP GND EO. SELF’) GOTORETG 


100 W = WHV1(WING(KL,10),WING(KL,12) ,WING(IJ,6),WING(IJ,8)) 
* - WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,9)) 
* - WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6),-WING(IJ,8)) 
* + WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6),-WING(IJ,9)) 
W = W*.25/PI 
RETURN 
110 W = WHV2(WING(KL,11),WING(KL, 12) ,WING(IJ,6),WING(IJ,8)) 
* ~ WHV2(WING(KL, 11) ,WING(KL,12) ,WING(IJ,6),WING(IJ,9)) 
* - WHV2(WING(KL,11) ,WING(KL, 12) ,WING(IJ,6),-WING(IJ,8)) 
* + WHV2(WING(KL,11),WING(KL,12) ,WING(IJ,6),-WING(IJ,9)) 
120 W = W*.25/PI 
RETURN 
END 


WAXHAW ek a wesw ak Kk ok sek ak ae ok ek wh ak ok aE aE ah ak ak ok ak ab Ok oe ae ok ak ok ok ak ok oe ok ab ok ok ok ok ak ak ah ak ak ak Ok ee ak OK oe aE ah ok ak ae ae Ok OE IK Ok 

* * 

* FUNCTION WHV1(X1,Y1,X2,Y2) 

x* eae 

bk bbe k bok kb obok bbb ik bok bb bb bk bbb b bok ak shot ak ab ab ak ob sh ak ok aba ob ab ok ak ok oh ok ob ob ak ah ok ak oh Ok ak a oe ok ah OF ok YK oe ok 
FUNCTION WHV1(X1,Y1,X2,Y2) 


IF (ABS(X1-X2).LT..0001) GOTO 100 


WHV1 = (1.0+SQRT((X1-X2)**2 + (Y1-Y2)**2) /CKHEK2) Gy ee 
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RETURN 
HOO” Wives 1 .0/(Y1-Y2) 
RETURN 


END 


Keak kkk Kk kk ak ek kak ak ak kk ak ak ak ak sek K HK Me Hk He kk ak kk kk kk kak Kak kk ak kkk kK KKH eK HK 
* * 
* FUNCTION WHV2(X1,Y1,X2,Y2) 
x* * 
KK KKK KKK Kk kkk Kk kk KKK Kk kK kkk kk Kk ok ke sek ak ak ak ak ak ak ak Mk a ak ae ak ak Ok ak ok ak sk ak a ak ak ak Uk a ok Mk 
FUNCTION WHV2(X1,Y1,X2,Y2) 
MENCABS(X1-X2 em. 000M) GOTO 100 
WHV2 = (1.0+(X1-X2) /SQRT((X1-X2)**2 + (Y1L-Y2)**2))/(Y1-Y2) 
RETURN 
100 WHV2 = 1.0/(Y1-Y2) 
RETURN 


END 
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MULT 1 


KKKKKKKHKKKKKKKKHKKHK KK HK KKK KKK HKKKHK KKK KKK KH AK KK KKK KKK KEKE KK KEK KE KEN EK KKK EKKEK 


a 


* 


SUBROUTINE MULT1 


SUBROUTINE TO MULTIPLY THE SPECIFIED LOADING VECTOR, 


(i 


BY THE COEFFICIENT MATRIX, AAC ), TO GET THE WING SHAPE 


VECTORVCL. 13) 2 


VARIABLE NAMES ARE CONTAINED IN THE MAIN PROGRAM. 
ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON. 


THE ONLY VARIABLE CHANGED IS WING(I,19) 


* 


* 


* 


KRKRERKEKEK KKK KKKEKKKEKERK EEK REE ERE KEE KK KK KKK KARE KEK EK KEKE KEKE KK KKEKKEEKK 


100 


SUBROUTINE MULTL 


REAL LAMDA 


CHARACTER*3 ELLIP 


COMMON PI,B,WING(101,23), SECTN(11,10) , AR, LAMDA, DELTA,NX,NY, 


* ALPHA,CLDSRD,ELLIP,XAC,CLA 


COMMON/COF/A(101,101), NEQNS 


COMMON /CAM/AA(101,101), AS(101,101) 


DO 110 I = 1,NEQNS 
WING(I,19) = 0.0 


DO 100 J = 1,NEQNS 


WING(I,19) = WING(1I,19)+AA(I,J)*WING(J,17) 


CONTINUE 
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MULT 2 


KEKKKKKAEKEKKAKA KEK KEK EEK KKK KER AERA KEKE K EKER EKER AAA KEKE EKER ER EK KEK ERK 


<< * 
* SUBROUTINE MULT2 * 
* * 
* SUBROUTINE TO MULTIPLY THE FLAT PLATE LOAD VECTOR, (1,15) * 
* BY THE COEFFICIENT MATRIX, AS( ), TO GET THE DOWN WASH * 
* VECTOR Gl 22.0 * 
* * 
* VARIABLE NAMES ARE CONTAINED IN THE MAIN PROGRAM. * 
* ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON. * 
* THE ONLY VARIABLE CHANGED IS WING(I,22) * 


* * 
KEKKEKKEKEKEKRKEKKER KEKE EK KEK KE EKER EE KEK KEE KK EE KEKE EE KEKE EEK EKER EK EE EATERS YS 
SUBROUTINE MULT2 
REAL LAMDA 
CHARACTER*3 , ELLIP 
COMMON PI,B,WING(101,23), SECTN(11,10),AR, LAMDA, DELTA, NX,NY, 
* ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101), NEQNS 
COMMON/CAM/AA(101,101), AS(101,101) 
DO 110 I = 1,NEQNS 
WING(I,22) = 0.0 
DO 100 J = 1,NEQNS 
WING(1I,22) = WING(I,22)+AS(1,J)*WING(J,15) 


100 CONTINUE 


14° 


a 
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MU 3 


KKK KHH HK HK HHH KH HH HHH HK HHH KH HHH HK HHH HHH KWH HHH HH MH H HH WK RHA HH HHH HH He MK HHH 


* * 


* SUBROUTINE MULT3 * 
* * 
* SUBROUTINE TO MULTIPLY THE SPECIFIED LOADING VECTOR, (1,17) * 
* BY THE COEFFICIENT MATRIX, AS( ), TO GET THE DOWN WASH * 
* VECIOR .Cl-25) * 
* * 
* VARIABLE NAMES ARE CONTAINED IN THE MAIN PROGRAM. * 
* ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON. * 
* THE ONLY VARIABLE CHANGED IS WING(I, 23) * 


* x* 
KEKEKEKKE KK KEKE KK KKK HHH HKK HK ke 9k Hk ke Hk HK H ake vk ve Heke HK Hk Hk Hk ak ok Hk ek Kk Hk KK Hk kK kK Kk kK KKK 
SUBROUTINE MULT3 
REAL LAMDA 
CHARACTER*3 ELLIP 
COMMON PI,B,WING(101,23), SECTN(11,10),AR, LAMDA, DELTA,NX,NY, 
* ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101), NEQNS 
COMMON /CAM/AA(101,101), AS(101,101) 
DO 110 I = 1,NEQNS 
WING(I,23) = 0.0 
DO 100 J = 1,NEQNS 
WING(I,23) = WING(1,23)+AS(I,J)*WING(J,17) 


100 CONTINUE 
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GRAP85 


KEKKKERK KKK KKK KKK KKK KK KE KEKE KK KEE KEK KEKE EK KEKE KEK ER KKK KKK EK KKK EERE KEK EEK 


* 


* 


PROGRAM GRAP85 * 


THIS IS A PROGRAM WHICH GENERATES A MENU OF GRAPHIC OPTIONS * 


AND WRITES FILES TO DISK WHICH ARE USED BY GRAPHER, A SOFTWARE * 


PROGRAM FROM GOLDEN SOFTWARE, INC. BOX 281, GOLDEN, COLO., * 
80402. IF YOU DON’T HAVE A COPY OF THAT SOFTWARE, YOU MUST * 
USE SOME OTHER METHOD OF PRESENTING THE DATA. * 
* 

DATA FILES USED FOR GRAPHIC OUTPUT ARE: * 
30 = X,Y DATA FILE * 
32 = GRAPH FILE * 
* 

* 

ALL PARAMETERS ARE PASSED IN COMMON BLOCKS. * 
VARIABLES USED ARE: * 
* 

ELLIP : ‘YES’ OR ‘NO ', DEPENDING ON WHETHER THE SHAPE * 
WAS COMPUTED a 

TYPE : GHARACTER STRING "FLAT ‘ OR ’CAMBERED’ * 
ITYPE : INTEGER, 1 IF FLAT, 2 IF CAMBERED * 
TITLE : THE TITLE OF THE GRAPH, 6 POSSIBLE TITLES ARE USED * 


% 


ISELEC : INTEGER RESPONSE FROM THE GRAPHIC SCREEN 


XMIN : MIN VALUE ON HORIZONTAL AXIS FOR PLOT OF A SECTION * 
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YMIN 


HS TART 


VINC 


AXINC 


YINC 


PLE 


alia 


XLE 


XTE 


AXIS 


TITLE 


Par Le 


je gela 


MAX VALUE ON HORIZONTAL AXIS FOR PLOT OF A SECTION * 


MIN VALUE ON HORIZONTAL AXIS FOR PLOT OF A SPAN 


MAX VALUE ON HORIZONTAL AXIS FOR PLOT OF A SPAN 


MAX VALUE ON VERTICAL AXIS OF PLOT 


MIN VALUE ON VERTICAL AXIS OF PLOT 


A TEMPORARY VARIABLE USED TO HOLD A VALUE. 


AN OFFSET USED TO LOCATE THE CORRECT COLUMN IN THE 


WING MATRIX 

SPANWISE SECTION TO BE PLOTTED 

THE ESCAPE CHARACTER (27) 

POSITION ON THE GRAPH WHERE THE HORIZONAL AXIS 
STARTS 

VERTICLE AXIS INCREMENT IN UNITS 

HORIZONTAL AXIS INCREMENT WHEN PLOTTING CHORD 
HORIZONTAL AXIS INCREMENT WHEN PLOTTING SPAN 
PLOT LEADING EDGE POSITION IN INCHES 

PLOT TRAILING EDGE POSITION IN INCHES 

LEADING EDGE POSITION IN PLANFORM UNITS 
TRAILING EDGE POSITION IN PLANFORM UNITS 
CHARACTER ARRAY WHICH HOLDS THE TITLES FOR THE 
HORIZONTAL AXIS. 

CHARACTER ARRAY WHICH HOLDS THE TITLES FOR THE 
VERTICAL AXIS. 

CHARACTER ARRAY WHICH HOLDS PART OF THE GRAPH 
TIOLE. 


CHARACTER ARRAY WHICH HOLDS PART OF THE GRAPH 
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* 


pad 


* 


* 


COEF 


COEFF 


SUBROUTINES 
CDATES = 


DELAY , 


SDIG 


’ 


RNDUP , 


PETCE 

CHARACTER ARRAY WHICH HOLDS PART OF THE GRAPH 
led ba by 

REAL ARRAY WHICH HOLDS THE LIFT, DRAG OR MOMENT 


ASSOCIATED WITH A PARTICULAR GRAPH. 


CALLED ARE: 
FORTRAN, GENERAL DATA ENTRY 
FORTRAN, SLOW DOWN THE PROGRAM SO ERROR MESSAGES 


CAN BE SEEN BEFORE BEING ERASED. 


FORTRAN, USED TO FIND THE NUMBER OF SIGNIFICANT 
DIGITS IN THE AAIS LIMITS ANDSMAKE AXIS Ease 
MORE UNIFORM. 

FORTRAN, USED TO ROUND THE AXIS LIMITS TO THE NEXT 


WHOLE DIGIT. 


* 


KEKE EK EEK AK AEE ALE AERA KERR RE ER AREER KER ARE ERR EERE KEE RR REEL AREER RSE ERR EES 


PROGRAM GRAP85 


REAL LAMDA 


CHARACTER*3 


ELLIP 


CHARACTER*8 TYPE(2) 


CHARACTER*20 TITLE(6) 


CHARACTER*8 TITLE2(2,2) 


CHARACTER*1 


ESC 


14° 


CHARACTER*23 AXIS(2) 

CHARACTER*7 COEF(2,3) 

REAL COEFF(2,4) 

COMMON PI,B,WING(101,23), SECTN(11,10) ,AR, LAMDA, DELTA,NX,NY, 
* ALPHA, CLDSRD,ELLIP, XAC,CLA 

COMMON/COF/A(101,101), NEQNS 

COMMON /CAM/AA(101,101), AS(101,101) 

COMMON /FILS/LIN, IOUT, JOUT , KOUT, LOUT , MOUT , NOUT 


KRKEKKEKKKKEKKEEEKKKEKE KKK K EK KEKE KKK KKK KEKE KERR KERR RE KKK RK ERK EKER ERK ERK KK KEE 


* * 
* DEFINE THE CHARACTER STRINGS NEEDED FOR GRAPH TITLES * 
* 7 


KEKEKKKEKKKKKKKEKKKKEKKEK KEK KKK KKK ERE KEE KEKE KEE KK EK KEE KEE KEKE KE EE I 

TYPE(1)=' FLAT’ 

TYPE(2)=' CAMBERED' 

TITLE(1)= 'd(CL)/dy é' 
TITLE(2)= 'd(CD) /dy ' 
TITLE(3)= 'd(CMAC) /dy ' 
TITLE(4)= ‘TWIST IN DEGREES ' 
TITLE(5)= ‘DELTA CP ' 
TITLE(6)= ‘WING HEIGHT COORD Z ' 


Pithe2(1,l)= ", AOA = * 


BITLE2(1,2)= ° deg 


MITLE2(2,1j)=] ', CLi = ° 
TITLE2(2,2)= ' ' 


AXIS(1)= ‘'SPAN-WISE COORDINATE Y ' 
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AXIS(2)= ‘'CHORD-WISE COORDINATE X’ 


COEF(1,1)=‘CL 


COEF(1,2)='CD 


COEF(1,3)='CMAC 


COEE(Z 1) — cles 


COEF(2,2)='CD 


COEF(2,3)='CMAC 


¢ 


BREE KKKEKHKHKAKHHKAK KKK KK KAKA AKA KEKE KKK KANE KEK KHKHEAE KEK EKEKEKAKKE KEE EE KEKEKERKERES 


* 


* ASSIGN DEFAULT VALUES TO THE VARIABLES 


* 


x 


* 


wv 


KAKEAEAKEKKAKKAKEKEKK KEE RK KEE KEKE KK KEKHKKEK KEK KEKE EER EEE SK LE KE KK KK HK KH KHA KEKE KE KK 


IMOD = 1 
AR= 2.0 
LAMDA = 1.0 
ADELTA= 1.0 
NX= 4 
NY= 4 
ALPHA = 5.0 
ELLIP = ‘YES’ 
CLDSRD= 1.0 
ITYPE = 1 
ISECT = 1 
XAC = 0.0 


ESC = CHAR(27) 


PI = 2.0*ACOS(0.0) 
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FARIA HHH HHA HAHAHAHAHA HHH HHH HHH HH HHH EHH HER HEHE HEHEHE HEE DEE FEM FE HEE FE KEE EE HEE HEE HE 


* * 
* SET INPUT/OUTPUT VARIABLES AND OPEN THE NECESSARY FILES * 
* * 


FH KAKA KKH KHAKI EHH IE EE EE IEE IEE EE EE IE IE IE I IE DE IE IEE TE DE TE HE IE TE IE IE OE OE EK 


IIN = 5 

IOUT = 6 

JOUT = 10 

KOUT = 11 

LOUT = 12 

LDAT1 = 13 

LDAT2 = 14 
OPEN(UNIT = IIN) 


OPEN(UNIT = IOUT) 

OPEN(UNIT = JOUT, FILE=’WING.MAT’ ) 

OPEN(UNIT = KOUT, FILE='’SECTION.MAT’ ) 

OPEN(UNIT = LOUT, FILE='PARAMS.DAT’ ) 

OPEN(UNIT = LDAT1, FILE=’FLAT.DAT’ ) 

OPEN(UNIT = LDAT2, FILE=’CAMBER.DAT’ ) 
10 FORMAT (16/3(F6.2,/),2(13/),F6.2/A/F6.2/12/12/F6.2) 
20 FORMAT (/2(10X,F12.6/)/10X,F12.6) 


KKKKK KEKE KEKE KEKE KI KIRKE KKH K KEKE KIER ER EK ERK KERR EERE KEKE REESE ERE REEL ERE 


* * 
* FORMAT STATEMENTS ARE FOR THE INPUT TO THE GRAPHER PROGRAM. * 
* OTHER GRAPHICS PROGRAMS WILL REQUIRE DIFFERENT FORMAT STMTS. * 
* * 
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KAKKKKKAKKKERK ERE KEK AK ERA EAE KKK KKK KAKA AK KAKA NK KK KAN KAA AK AK KEK KEKE KK KEKE 


32 


Se 


34 


FORMAT ( 

* 11236'/ 

x eo 1 Oo 

* PL 

* '65 66 78 "NO" O'/ 

* '"NO" "SOLID" 1.500e-001 1'/ 

* 'YES" 41 1.000e-001 1 1'/ 

* '78 9.900e+028 9.900e+028 0.000e+000 "DEFAULT" 1.000e-001 1'/ 

* 'SOLID" 5 1.500e-001 0.000e+000 9.900e+029 100 2.000e+000 1’) 
FORMAT ( 

* '1236'/ 

Soi eon ay 

cot ploy 

* '65 66 78 "NO " 0'/ 

* ' NO" "SOLID" 1.500e-001 1'/ 

* "YES" 41 1.000e-001 1 1'/ 

* '78 9.900e+028 9.900e+028 0.000e+000 "DEFAULT" 1.000e-001 1'/ 
* SOLID" 5 1.500e-001 0.000e+000 9.900e+029 100 2.000e+000 1°) 
FORMAT ( 


* 'HORIZ! / 


% 


'1.750e+000 ',E10.3,’ 6.000e+000 88'/ 
* E10.3,’ ',E10.3,’ 9.900e+028 1 1'/ 

* 10.000e+000 ',E10.3,' 1.500e-001 1 1'/ 
+ ey 


* 1 2.750e+000 4.000e-001 9.900e+028 1.800e-001'/ 


2 


36 


5] 


38 


oo 


* '"DEFAULT" "DEFAULT" "',A,‘"') 
FORMAT ( 
SO yeARTS' / 
* '1.750e+000 1.000e+000 5.000e+000 89'/ 
* E10.3,' ',E10.3,’ 9.900e+028 1 1'/ 
* '2.700e+002 ',£10.3,' 1.500e-001 1 1'/ 
ee 11 ail / 
* '1 9.900e+028 0.000e+000 9.900e+028 1.800e-001'/ 
* '"DEFAULT" "DEFAULT" "',A,‘"') 
FORMAT ( 
* PLY 
* '"DEFAULT" 1! / 
* '2.000e+000 7.000e+000 0.000e+000 1.700e-001' / 
k O'/ 
* °0 13 "' A,’ WING"’/ 
cee 12 "CLi = *,F5.2,0"') 
FORMAT ( 
~ Pi! / 
* ' "DEFAULT" 1! / 
* '2,000e+000 7.000e+000 0.000e+000 1.700e-001'/ 
soe)? / 
oO 13 "* A,’ WING"’/ 
pe 34°" A F6 149A, F5.2,A'"*) 
FORMAT ( 
Zep” / 
* ' "DEFAULT" 1// 
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* '2,.000e+000 7.000e+000 0.000e+000 1.700e-001' / 
® O'S 
* 'O 26 "°. A,” WING, SECTION. sailor sane. 
* '] 26 "vertical scale exaggerated"’ ) 
40 FORMAT ( 
x TE 7 
oe DE A ee lee 
* '2.000e+000 7.000e+000 0.000e+000 1.700e-001'/ 
ee) 
xO 26%" A, WING) SHCRION: “les ) 
41 FORMAT ( 
* LE 
co SDE RAUIORYy 1827 
* E10.3,’ 6.000e+000 0.000e+000 1.000e-001'/ 
eee 
a UO Ra. Miss) 
42 FORMAT ( 
saa) D4 
i DE PAU ICI Sk 7 
* E10.3,' 6.000e+000 0.000e+000 1.000e-001'/ 
nay 
-AO yO ee) 
43 FORMAT ( 
come ded! 
k'DT/ 


* E10.3,’ 1.000e+000 ',E10.3,’ 6.000e+000 1 2 1.000e-001'/ 


‘lab 


&~4 


FEO ,3,°)1,000e+000 * ElOT3," 6,.000e+000 1 2 1.000e-001") 
50 PORMAT (Z2(F10.6 1X) > 
60 FORMAT (23F10.6) 


70 FORMAT (10F10.6) 


KK KKE HK HH HH HH HHH HH HHH HAH HHH HWW HH WH HW EW HE Ek I dk BE I dk k ak ak Dk ake IE IE IE ak ak ak Hk He MEK HK 


* * 
* THESE FORMAT STATEMENTS ARE FOR THE FLASHUP WINDOWS MENUS * 
* * 


Kk kkk kk kk kk kk kk ok kk kk ek Hk KH HK ak He Ik ok ok ak ak Hk Dk Ok ok ok ak Hk ak ak ak Ok OK kk Hk ak Hk ak ak Ik KE ME ek kk HK KH HK HK 
82 FORMAT (' ',A,A,13) 
84 FORMAT (' ',A,A,A) 


KRAKEKAEKEK AKER KEKE KEKE KKK EER EEE REE ER EERE EERE ERK EEK EEK HEH HK KEKE EK EK 


* * 
* READ THE LAST SET OF PARAMETERS, IF UNCHANGED, READ WING AND * 
* SECTN DATA , IMOD IS THE FLAG TO TELL IF THE PARAMETERS HAVE * 
* BEEN CHANGED. 1 = CHANGE, 0 = NO CHANGE. * 
* * 


HAHAHAHA HHH HH RH HHH HHA HAHAH HHH HHH HHH HAA AAA HH HHH HHH WIE TE DE BEE ak ak ak ak ak ak ake ak ae 
READ (LOUT,10) IMOD,AR,LAMDA,ADELTA,NX,NY,ALPHA, 
ZOELLIP,CLDSRD,ITYPE, ISECT,XAC 


DELTA 


ADELTA*PI/180. 


NEQNS = NX*NY 


B = AR*(1.0+LAMDA) /2.0 
IF (IMOD.EQ.0) THEN 
READ (JOUT,60) ((WING(I,J) ,J=1,23),1=1,NEQNS) 


Rear cCKOUr, 70) “CCobCIN(!, J) J=1,10),1=1r,NY) 
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READ (LDAT1,20) (COEFF(1,J),J=1,3) 
READ (LDAT2,20) (COEFF(2,J),J=1,3) 
ELSE 
WRITE (IOUT,*) ' CHANGES WERE MADE TO THE PLANFORM. GO BACK’ 
WRITE (IOUT,*) ‘ AND RUN THE LOAD/SHAPE PROGRAM AGAIN. ' 


CALL DELAY (6) 


GO TO 1000 
ENDIF 
COEFF(1,4) = ALPHA 


COEFF(2,4) = CLDSRD 
CLOSE (JOUT) 
CLOSE (KOUT) 
CLOSE (LDAT1) 
CLOSE (LDATZ2) 


KKKRKKEKKRREEERERERRRERKREEERKEERERERRRERRRREREREREREREREERESRKER ST Oe 


* * 
* SET MINIMUM VALUE OF X. THIS X IS THE CHORDWISE COORDINATE * 
* ¥* 


KHAKHKK KKK KKK KK KK KK AK K KKK KKK KKK KK KK KKK EK KEK KE KKK KEE KKK KEK EK EE KKK KEE KK KEK 


XMIN = 0.0 


KEKKAKEK AK KKK KKK KEKE LEE KEK KK KKK EK KK KEK KEK HK KKK KEK KE KKK KEE KKK KE KK EKG KEKE 


mo a 
* COMPUTE THE CORNER OF THE WING, TO GET MAX VALUES FOR THE AXIS * 
* THIS IS NEEDED IN CASE THE WING IS SWEPT OR HIGHLY STAPERED: * 
* * 


KKKKKKKK KKK KKK KK KKH KE KEKE KKK KEK EK EEK KKK EK HK EK KEKE KEK KKK EK KEE KK KE RL EE KKK EK 


TSG 


XT = (B*TAN(DELTA) /2.0) + LAMDA 


KREKEKEKAK KEKE AKA KE KEKAK AK KEKE KK KEKE REE KEKE KER ERR K KK KEKE KKEKKEKERERER ERR EE ERE 


* * 
* NOW CONVERT SPAN AND CHORD AXIS LIMITS INTO "ROUNDED" UNITS * 
* Vitaeeoane TIME, COMPUTE THE INCREMENTS IN TIC MARKS FOR THE * 
* GRAPHS. * 
* * 


HREKEEAEKK KEKE KEKE KKK KEKE KEKE KKK KEK ERK KEKE RERERR EE RKERERRERRRRRRRR RRR ERE 


ieGxaecT.1,) THEN 
CALL RNDUP(XT, XMAX) 


ELSE 


ENDIF 
YMIN = 0.0 
CALL RNDUP(B/2.0, YMAX) 


XINC 


(XMAX -XMIN) /5.0 


YINC 


(YMAX-YMIN) /5.0 


HK KAKEEER ERE EEK KEKE KEKE EEE ERE EKER EK EERE REE ERE EEE EEE EEE EEE OE OIE IF 


* * 
* CLEAR THE SCREEN AND WRITE THE GRAPHIC MENU. * 
* * 


KAEKREKAKEKA KAA EKER KEKE KEE KK EEE EO i bbb bb bk bk bt ot kkk bt 
110 WRITE (IOUT,*) ESC,’[2J' 
WRITE (IOUT,*) ESC,’[00;00H’ 


WRITE CLOUT *) " C=ALL/” 


hoy 


WRITE (IOUT,*) ' W=GRAPH/’ 
WRITE (IOUT,84) ESC,’[08;69H’ , TYPE(ITYPE) 
WRITE (IOUT,82) ESC,'[09;69H’ , ISECT 


KEKKKKKKEK KK KEKE KKKEKKKE KE KK KEK KKK KE KKKE KEKE K KKK KEE KK KEK KK KK KK KK KKK KKEK KEK 


* * 
* READ THE OPTION AND SEND TO THE APPROPRIATE AREA. * 
* OPTIONS ARE: cs 
* 1, SETS THE TYPE OF WING, FLAT OR ELLIPTIC. * 
* 2, SETS THE WING SECTION TO BE PLOTTED * 
* 3, RETURNS TO THE MAIN MENU * 
* 4 - 7 PLOT THE LIFT, DRAG, MOMENT OR TWIST VS SPAN * 
* 8 PLOT DELTA CP VERSUS CHORD * 
* 9 PLOTS THE WING SHAPE VERSUS CHORD, WHICH GIVES * 
* AN ELLIPTIC LIFT DISTRIBUTION. * 
* * 


KRKEKKKKK KEKE KE KKKEK KEK KKK KK KKH KKK KK HHH KK HHH HH KK KK KK KK KAA KEK KKK AN KKNEN SE 

READ (IIN,*) ISELEC 

WRITE (IOUT,*) ESC,'[2J' 

WRITE (IOUT,*) ESC,’ [00;00H' 
KEKE KKEK KE KEKE KKK HH EK EK HHH KH HHH HH HHH HH HHH HH HHH MHHKAH RH KKH HHKHKKKHK KKK KKK KKK AK 
* + 
* THIS IS OPTION 1, SET TYPE OF WING. * 
* x 
KKK KKK HH KKK HH KKH HHH HH KKH H HHH HHH HH MHKK HHH HHH HH HHH HH HH HK HK HH KK KKK KK KKK KKK K 


IF (ISELEC.EQ.1) THEN 


Tee LLIF SEOs NOe) Shar 
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WRITE (I0OUT,*) ‘' ELLIPTIC DATA WAS NOT COMPUTED’ 
CALL DELAY (6) 
GO TO 110 

ELSE 

ENDIF 


weecliyPEe EO: 1) THEN 


DO Oi 2 ae 
ELSE 
TYEE <= — 1 
ENDIF 
GO TO 110 
ELSE 
ENDIF 


HKKKKKKKKKKHKK KEK WK KKH KKK KKH WK KKH KKK HK KHNW HH KK HW WW HH HK MWe WH HK KW MW HH Hk KKH HK HK Hk KM" HK HK 


* * 
* TvioelosOrrl@N 2, SET SECTION NUMBER * 
* * 


HHKKK HHH HHH HAHAH HHH HH HHA KAHL HHH HAH HHH KAM H HH HHH HHH HK HH HHH HH HE HK 9K 
DEP OrSELEC. EQ. 2) THEN 


CAE CUTTGas( LyNY , 1SECT, 


* ‘WING SECTION NUMBER ’ , IMOD) 
IMOD = 0 
GO TO W110 
ELSE 
ENDIF 


KREKEKKKEKEKKE KKK KE KEK KK KKK HHH HHH HH HHH HHH HH HH DE EDK DE TH IK DK DK OK Ok ak Ok ak Ok ok ok ak ok ak ok ok ok ok ae AK ok ok ok OF 


5 9 


* * 


* THIS IS OPTION 3, RETURN TO MAIN PROGRAM. BEFORE LEAVING, WRITE * 
* THE FLASHUP WINDOW FOR THE MAIN MENU. * 
* * 


KKEKKKKKKKHKKRKKKEKRK KK KK KKK KKK KKK Kwek Kk KK Kk Kk HK ok ok kkk kk Kk Kk Kk kkk Kk Kk KK KH HK 
IF (ISELEC.EQ.3) THEN 
GO TO 1000 
ELSE 
ENDIF 
IF ((ISELEC.LT.1).OR.(ISELEC.GT.9)) THEN 
WRITE (IOUT,*) ‘' THAT IS NOT A VALID SELECTION’ 
CALL DELAY(6) 
G0 107110 
ELSE 
ENDIF 


HHKKKKKKKKKE KEKE KEKE EKER KKK KEKE KEKE KER EK ERE KER ERE KEKE KEK EK EK EKER EEK EK EKER 


x * 
* THESE ARE THE LIFT, DRAG , MOMENT AND TWIST VS "SPAN PLEOTS= * 
* * 


KEKEKKEREKEKKEKEK KEKE EK EK KE KEKE REE KEKE EKER EE EKER EE KKK KEE KEKE KEK KEE RK KEK 


LP aC SELEC CQ.) reO, 10.55 10 


FKRKKKEKEKE KEK KKK KEKE KEK KEE KE KEKE KEK KEK EK RK KEKE KE KEK KKK KK KK KEK KEKE KKK KEK EK EEK RERER KEK 


* * 
* OFFSET IS ESTABLISHED, IT IS AN EASY WAY TO USE THE StsEe eas * 
* TO FIND THE CORRECT COLUMN OF THE SECTN MATRIX. * 
* SINCE THERE IS NO TWIST ON A FLAT WING, SET THE WING TYPE TO * 


La! 


k= 


* CAMBERED IF THE TWIST PLOT IS SELECTED. * 
* * 
FEE EAE IEE HEE TE TEE IE DE EE DE TE EE DE HEE DE EE EEE ETE HEE HEH EHR HEH HEA HEA HEE ETE BEE HE DE TEE DE ETE HE WIE KE KE 
iP CESELEC.EQ.7) @MIYPE = 2 
Pie clul wer EQ. 1) THEN 


IOFFS -2 


ELSE 


I 
Ww 


IOFFS 


ENDIF 


HHH HH KK KK HHH HHH KK HK HHH HE HK HHH HK HH HK HH IE HH HM HHH HE EE He ME IE He IE He EE KE EEE KE Fe aE IEE ee Ee EK 


* * 
* SET UP THE DATA FILE FOR THE GRAPHIC PLOT. * 
* * 


- 


HHH KKKAEK KEK KK KKK KEKE KKK KK KK KKK HK KK KKK KKK KK KKK KK KK KK KK KEKE KEE EEE LE EEN 


sf 


J=1 
VMAX = -100.0 
VMIN = 100.0 


OPEN(UNIT = 30,FILE ='P1.DAT’) 
IJ = J 
BOeS00 L = 1,NY 

IJ = (I-1)*NX + J 

VERT = SECIN(1, IOFFS+ISELEC) 


IF (VMAX.LT.VERT) VMAX = VERT 


IF (VMIN.GT.VERT) VMIN = VERT 
WRITE (30,50) WING(IJ,5), VERT 


300 CONTINUE 
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CLOSES) 


KKEKKKEKKKKKKKEKKEKE KERR EER EE REE KER KEKE EK KRE KEE KKK KKK KR KKKKKRKK RRR RRR ERRRRRE 


* * 
* SET UP THE P1.GRF FILE FOR GRAPHER. IT HAS ALESTHEeGRer * 
* FORMATS IN IT. * 
* * 


KKKKE KEKE KKK K KEKE EK KEK KKK ER KEK KEE KER EEK KERR REE K KKK ER EK EK KK EER KE KEK KEK ER EEE 


IF (VMIN.GT.0.0) VMIN 070 


IF (VMAX.LT.0.0) VMAX = 0.0 
CALL RNDUP (VMAX, VMAX) 
CALL RNDUP(VMIN, VMIN) 
CALL SDIG(VMIN, VMIN, VMAX, VMAX, IDIG) 
VINC = (VMAX-VMIN)/5.0 
HSTART =(VMAX-6*VMIN) /(VMAX-VMIN) 
OPEN(UNIT = 32,FILE ='P1.GRF’) 
WRITE (32,32) 
WRITE (32,34) HSTART, YMIN, YMAX, YINC, AXIS(1) 
WRITE (32,36) VMIN, VMAX, VINC, IDIG, TITLE(ISELEC-3) 
IF (ISELEC.EQ.7) THEN 
WRITE (32,37) TYPE(ITYPE) , COEFF(2,4) 
ELSE 
WRITE (32,38) TYPE(ITYPE) , COEF(ITYPE, ISELEC-3), 
*  COEFF(ITYPE,ISELEG-3), TITLE2(1TYPE, lj COB Fr Clement 
*  TITLE2(ITYPE, 2) 


ENDIF 


CLOSE 72) 
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GO TO 1000 


FOKKER AAAI HHH KKK HEH K KKK HK HHH HK HK HH HEHEHE HEHE FEE HE EE HE FEE HEHEHE FEE HAH FE KH HHH HEHEHE KH 


* * 
* THESE ARE THE DELTA CP PLOTS * 
* * 


SKAK AKA AA AAA KHAN NAH HAA AKAIKE AE HEE HEHE EE EVE DEE HE EEE DEDEDE TE HEHE HEH HEM PEW He 
310 CONTINUE 


IF (ISELEC.EQ.9) GO TO 400 


HKKKKKKK KKK KKK KK KEK KEK EK KE KKKEE KEK ERK ERK KE KKK ERK EKER EKER KEKE RK KEK EK KERK EK 


* * 
* OFFSET IS SET AGAIN. x 
* * 


KEKKKKKKKKK KKK KKK KKK KK KK KEK KK KEKE KKK KE KEKE RK KEKE ERK KKKKK KEKE EK KKK KKK EEE EK 


IF (ITYPE.EQ.1) THEN 


IOFFS = 8 
ELSE 

IOFFS = 10 
ENDIF 


BHKKKHKA AA KKK HA KKK KEK AK HHH KKK KK KEK EK KEKE KER KE KEK AE REESE KEKE KEK KK EN 


* * 
* SET UP THE P1.DAT FILE FOR GRAPHER iS 
* 


* 


BHKKKAK KER EEA KKK HHH KKK KKK KEKK IKK KIKI EK IE HHH EEITEEE IEEE EEE TE E 


I = ISECT 
VMAX = -100.0 
VMIN = 100.0 
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OPEN(UNIT = 30,FILE =’P1.DAT’) 
IJ = (I-1)*NX + 1 


KRKKKBKEKKK KA KK NK KK KKK KK KKK KKK KK KKKHK KEE K KK KEK KE KKE KKK KEK KAN K EEK KKK EK EEEK 


* * 
* ALE IS THE VALUE OF THE LEADING EDGE. SINCE THE DEE Ser * 
* FOR THE CAMBERED WING IS ZERO AT THE LEADING AND TRAILING * 
* EDGE, IT IS POSSIBLE TO EXTEND THE DATA TO THOSE (Grit * 
* * 


WHEKKKKKEK KEKE KKEKEK KBE KBE EK EEK KEKE KKK KEK KEKE KEK KEK KEK KKK KK HK KKK KK HKA KKK KRK 
XLE = WING(IJ,4)-WING(IJ,13)/2. 
IF(ITYPE.EQ.2) THEN 
WRITE (30,50) XLE, 0.0 
ELSE 
ENDIF 


KHEKEKKK KE KEK KEKE KK KK KK KKK HKKKHK KKK WK HK HH WKF WK HK KWH HH WW HHH WH HH HWW HK WK HK HHH HH HK 


* * 
* GET THE MAX AND MIN VALUES OF THE PARAMETER BEING PLOTTED * 
* SOQ ALL GRAPHS WITH THAT PARAMETER WILL HAVE THE SAME VERTICAL * 
* SCALE. * 
* * 


KRKKEKEKK KAN KEK KE KKK KEKE KE HEH HH HK HK KKH HH HH EHH KH KKK HH HH HH HH HHH HRW KKK ke Hk ake WK 
DO 315 IJ = 1,NEQNS 


VERT = WING( IJ, IOFFS+ISELEC) 


IF (VMAX.LT.VERT) VMAX = VERT 


Ver 


IF (VMIN.GT.VERT) VMIN 


315 CONTINUE 
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KAKA KKK KKK KEK AKER KEK KE KEK KKK EK EK AE KARE KER ER ERE KEKE RE ERE RRR ER ER KEKE ERE EERE 


* * 
* SET UP THE PL.DAT FILE FOR GRAPHER * 
* * 


KKK KRE KKK KKK EK BRK KK RK KEKE KEKE AK KEAKKE KEK AK K KARE KAR KK KEKE AK EKA KEK KE AEE EK 
DO 320 J = 1,NX 
IJ = (I-1)*NX 4+J 
VERT = WING(IJ, IOFFS+ISELEC) 
WRITE (30,50) WING(IJ,11), VERT 
320 CONTINUE 
IJ = I*NX 


KRKEEKAEKKEKABEKEKEKADLEABER RRR KE ER ERE RRR ERR RRR REE RRR ERERRRERKRRRERERER ESR 


* * 
* XLE IS COMPUTED FOR THE REASONS STATED ABOVE. * 
* * 


HKREKKKKKKEEKE KEKE KERR KEK ERE BREE ERK KEK RRR REE KKK KK REAR KEKE ES 
XTE = WING(IJ,4)+WING(IJ,13)/2. 
IF(ITYPE.EQ.2) THEN 
WRITE (30,50) XTE, 0.0 
ELSE 
ENDIF 


GLOSE (30) 


KHKKKEKKEK KEK ERE KKK KER EEK ERR EKER EKER ERK KEKE RE BRK EEK RE RE ER ERE RRR ER REE 


* * 
* COMPUTE AND STORE THE VALUES FOR THE P1.GRF FILE USED BY GRAPHER* 
* * 
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KHHK WHR KHANH AKA AKKAKAK KKK AK KKH KK HHH HH KH HK KH HH HH HH HW KWH HHH HH HHH HK HH HH Hk Hk HK Hk Hk Hk HH 


IF (VMIN.GT.0.0) VMIN 0.0 


il 


IF (VMAX.LT.0.0) VMAX 0.0 


CALL RNDUP (VMAX, VMAX) 

CALL RNDUP(VMIN, VMIN) 

CALL SDIG(VMIN, VMIN, VMAX, VMAX, IDIG) 
VINC = (VMAX-VMIN)/5.0 

HSTART = (VMAX-6*VMIN) /(VMAX-VMIN) 


PE 


(6.*XLE/(XMAX-XMIN) )+1.75 


PAGES 


(6. *XTE/(XMAX-XMIN) )+1.75 

OPEN(UNIT = 32,FILE ='P1.GRF’) 

WRITE (32,33) 

WRITE (32,34) HSTART,XMIN, XMAX, XINC,AXIS(2) 
WRITE (32,36) VMIN, VMAX, VINC, IDIG, TITLE(ISELEC-3) 
WRITE (32,40) TYPE(ITYPE) , ISECT 

WRITE (32,41) PLE 

WRITE (32,42) PTE 

WRITE (32,43) PLE,PLE, PTE, PTE 

CLOSE (32) 

GO TO 1000 


KAKEKAK KK KKK KH HH KKH KKH HHH HHH HHH HHH HH KH H HH HEH HH HHH HH HH AH HH HHH HK HHA KH HK KKK KKK 


* * 
* THIS. £S THE SSHAPESELOT. * 
* * 


KKK KKK EKA KKK HK KEKE HK EKER KEKE ERK HEHEHE HHH HHH HHH HEHEHE HHH HE HHH HAW HHH HHH BE I 


400 IF (ELLIP.EQ.'NO ') THEN 


Life 


WRITE (IOUT,*) ‘ ELLIPTIC SHAPE DATA WAS NOT COMPUTED‘ 
CALL DELAY(6) 
GO TO 110 

ELSE 


ENDIF 


FR KKK RK KKK KAKA KKK KEKE KEK AHHH AK KK HK KKK KEK KKH HHH HK KKK HHH KKK KKK KKK IK KH 


* 


* 


pEE UP THE PL.DAT FILE FOR THE WING SHAPE. WE KNOW THAT THE 
LEADING EDGE STARTS AT ZERO HEIGHT, SO THAT STARTING VALUE IS 
ADDED TO THE DATA, IN ORDER TO GET THE MAXIMUM INFORMATION 
INCLUDED IN THE GRAPH 

XLE AND XTE ARE ALSO NEEDED SO THAT VERTICAL DOTTED LINES 


SHOWING THEIR POSITION CAN BE ADDED TO THE GRAPH. 


* 


* 


* 


* 


KAKAKAKEEKAKAKK AER AEA KKK KKK KK AK EKER KKK KEK KKK KARR ERK EKER KEK RK KEK KKK KKK 


t= ITSECT 
VMAX = -100.0 
VMIN = 100.0 


OPEN(UNIT = 30,FILE ='P1.DAT’) 
i Jee (1-1)*NX + 1 

XLE = WING(IJ,4)-WING(1IJ,13)/2. 
WRITE (30,50) XLE, 0.0 

IJ = I*NX 

XTE = WING(IJ,4)+WING(IJ,13)/2. 
DO 405 IJ = 1,NEQNS 


VERT = WING(IJ,20) 
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IF (VMAX: LT.VERT) VMAX "= "VERT 
IF (VMIN.GT.VERT) VMIN = VERT 
405 CONTINUE 
DO 410 J = 1,NX 
IJ = (1-1)*NX+J 
VERT = WING(IJ, 20) 
WRITE (30,50) WING(IJ,4)+WING(IJ,13)/2.0, VERT 


410 CONTINUE 


CLOSE (30) 
IF (VMIN.GT.0.0) VMIN = 0.0 
IF (VMAX.LT.0.0) VMAX = 0.0 


CALL RNDUP (VMAX, VMAX) 

CALL RNDUP(VMIN, VMIN) 

CALL SDIG(VMIN, VMIN, VMAX, VMAX, IDIG) 
VINC = (VMAX-VMIN)/5.0 

HSTART = (VMAX-6*VMIN) /(VMAX-VMIN) 


PLE 


(6.*XLE/(XMAX-XMIN) )+1.75 


PTE = (6.*XTE/(XMAX-XMIN) )+1.75 

OPEN(UNIT = 32,FILE ='P1.GRF’) 

WRITE (32,33) 

WRITE (32,34) HSTART,XMIN, XMAX, XINC,AXIS(2) 

WRITE (32,36) VMIN, VMAX, VINC, IDIG, TITLE(ISELEC-3) 
WRITE (32,39) TYPE(2) ,ISECT 

WRITE (32,41) PLE 


WRITE “(32542 PTE 


WRITES (32,43) PLE, PLE,PIE, PIE 
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CLOSE (32) 
1000 REWIND (LOUT) 


KXKAKAKE KEKE KEKE KEKE EERE KERR EKER EEE EER KEE REN ERE RRR ERE EER ERE ER ERE EK 


* * 
* RE-WRITE THE PARAMETER DATA FILE, SINCE IT KEEPS TRACK OF THE * 
* LAST WING TYPE AND SECTION PLOTTED. * 
* * 


KKK KKKKKK KERR K KK KKK KK HE KKK KKK EHH KEE HK KEKE KEKE KEE HEH KEKE KKK HK HK KEE 
WRITE (LOUT,10) IMOD,AR,LAMDA,ADELTA,NX,NY,ALPHA, 
*ELLIP , CLDSRD, ITYPE, ISECT , XAC 
CLOSE (LOUT) 


HAKEKKEKKKAAKKAERK AERA AERA AEA ERR ERE RRR ERE RERER ERR RRR RRR ERE RRR RERERRERER 


* * 
* BEFORE LEAVING THE PROGRAM, RE-WRITE THE MAIN PROGRAM MENU. * 
* * 


KAKKEEKERKEKEKEKE RK KKK KE RRR KKK R EERE REE ERK RRR REE RRR REE RRR ERR RRR RRRRERRE 
WRITE (LOUT,*) “ESC, ' [20” 
WRITE (IOUT,*) ESC,'[00;00H’ 
Were (rOur,*) * G=ALL/' 
WRITE (IOUT,*) ' W=MAIN/’ 
END 
SINCLUDE: ‘CDIT85.FOR’ 
SINCLUDE: ' DELAY. FOR’ 


SINCLUDE: ‘RNDUP. FOR’ 


SINCLUDE: 'SDIG. FOR’ 
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DELAY 


KHEKKEKKKEKKKEKERRERERE KER EERE KKK KEK KEK KKK KK KKK KKK EKKEKEKRRERKERRRKRKREK EERE 


* * 
* SUBROUTINE DELAY(CINTERV) * 
* * 
* CALLED TO SLOW THE SCREEN DOWN AFTER PRINTING AN ERROR MESSAGE. * 
* VARIABLES USED ARE: * 
* INTERV : NUMBER OF SECONDS FOR DELAY MUST BE LESS THAN 59, * 
* INEUS, * 
* IHR : TIME OF DAY IN HOURS * 
* IMIN ¢- SLIME OPS DAY SUN MIN * 
* ISECL : TIME OF DAY IN SECONDS * 
* ISEC2 ;: TIME OF DAY IN SECONDS * 
* ITHUND : TIME OF DAY IN HUNDREDTHS OF SECONDS. * 
* * 


KREKEKEKEKKEKREKERER EERE RRR EKER ERR REE RRR RRR EKER ERE RRERE ER ERERER 
SUBROUTINE DELAY(INTERV) 
INTEGER*2 IHR, IMIN, ISEC1, ISEC2, IHUND, INTER 
INTER = INTERV 
CALL GETTIM(IHR, IMIN, ISEC1, IHUND) 
100 CALL GETTIM(IHR, IMIN, ISEC2, IHUND) 
IF’ ((ISEG2-ISEG). LI SER) GODOM OO 
RETURN 


END 
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CDIT85 


KKK KEKK HAHA KKK KKK AAR KKK KKK AK KK KKK KKK KK EK KK KEK KE KKK KKK EERE ERK KKK KE KKK 


* * 
* SUBROUTINE CDIT85(LOW, HIGH, VALUE, NAME, IMOD) * 
* INTEGER VERSION OF CDAT85 * 
* HANDLE INPUT OF PARAMETER DATA FOR THE PROGRAM. INCLUDES SOME * 
* ERROR CHECKING. * 
* * 
* VARIABLES USED ARE: os 
* LOW : MINIMUM VALUE THAT WILL BE ACCEPTED, INPUT * 
* HIGH : MAXIMUM VALUE THAT WILL BE ACCEPTED, INPUT * 
* VALUE : VALUE READ FROM SCREEN, OUTPUT * 
* NAME Se CREE NPEROMPI STRING, INPUT * 
* IMOD : FLAG USED IN ANOTHER PROGRAM, SET TO 1 ANYTIME * 
* CDAT85 IS CALLED. , OUTPUT 


* ESC ; ESCAPE CHARACTER, (27) es 
* * 
KAKA A AHA HHA HHH A HHA HK KH A HHA IAA LAI A HEH AK HEATER HE ETE HEE TEEN TE EEE TER Be 

SUBROUTINE CDIT85( LOW, HIGH, VALUE , NAME, IMOD) 

COMMON/FILS/TIIN , LOUT , JOUT , KOUT , LOUT , MOUT , NOUT 

INTEGER LOW,HIGH, VALUE 

CHARACTER*35 NAME 

CHARACTER*1 ESC 

ESC = CHAR(27) 


10 FORMAT(’ ENTER THE ’ ,A35) 


ay 


20 FORMAT(’ THE MINIMUM VALUE FOR ’ ,A35,/’ IS ‘ ,13) 


30 FORMAT(’ THE MAXIMUM VALUE FOR ‘ ,A35,/' IS ’ ,13) 


KREKEEKKEEKKAKKKKER KEE KER KKK EK ERK EKRE RRR KEK EER REKKERKRERERERKEERKKEREREKKERE 


* * 
* CLEAR SCREEN, POSITION CURSOR, AND * 
* SET IMOD FLAG AND WRITE PROMPT TO SCREEN? * 
* * 


KEREKEKKERERKKEKKRRKRRRKERERREKRKKRERKERERERKKKRKKEXKREKKKKKKKKKKEKRKKEKAKR OH 
100 WRITE (IOUT,*) ESC,'[2J° 

WRITE (IOUT,*) ESC,'[10;0H’ 

IMOD = 1 

WRITE(IOUT,10) NAME 


KHEKKARAKAEEK EAE KK AKE KK KAK AK KK EKER KR KEE RE ERK RE KREK RE ERK KEK RK RE ERK RRR KEKE KER KERESR 


* 20S 
* READ SCREEN VALUE. ON ERROR OR END OF DATA SET, GO BACK TO * 
* PRONE T. * 
* * 


KKK AKKKEKRERR KEK AKER EKA KKK KKK ERK EK ERK ERK RERE ERK KE KEE RE RR EEK KEKE ER KERN E KES 


READ(IIN,*,ERR = 100, END = 100) VALUE 


KREKAKEKEKAKEEK ERK KKK KEKE RKEKER EKER RE RK ERK ER RRR KERR KEKE R RR KRE REE RREREERERERE 


* * 
* IF VALUES ARE OUT OF RANGE , PRINT ERROR MSG AND RETURN TO * 
* PROMPT. * 
* ¥ 


HKKEAKKKEKKKKKHK ERK KE KK EKA K AK KEK KEKE KEKE KEK ERE RK KEK EERE EKER KER RER ERR EER EN ERE 


IF (VALUE.LT.LOW) THEN 


WRITE (IOUT,20) NAME, LOW 
CALL DELAY (6) 
GO TO 100 

ELSE 


IF(VALUE.GT.HIGH) THEN 


WRITE (IOUT,30) NAME,HIGH 


CALL DELAY(6) 
GO TO 100 
ELSE 
END IF 
END IF 
RETURN 


END 
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CDAT85 


KXKKKKKKEKKE KEKE KEK EKER KEK KEE KEKE KKK KEKE KEKE RE EEK KERR KEK KEKE EKER KEK KE KEKE KE KKK ER 


* * 
* SUBROUTINE CDAT85(LOW,HIGH, VALUE, NAME, IMOD) * 
* REAL VERSION * 
* HANDLE INPUT OF PARAMETER DATA FOR THE PROGRAM. INCLUDES SOME * 
* ERROR CHECKING. * 
* * 
* VARIABLES USED ARE: * 
* LOW : MINIMUM VALUE THAT WILL BE ACCEPTED, INPUT * 
* HIGH : MAXIMUM VALUE THAT WILL BE ACCEPTED, INPUT * 
* VALUE : VALUE READ FROM SCREEN, OUTPUT * 
* NAME : SCREEN PROMPT STRING, INPUT * 
* IMOD : FLAG USED IN ANOTHER PROGRAM, SET TO 1 ANYTIME * 
* CDAT85 IS CALLED., OUTPUT * 
* ESC : ESCAPE CHARACTER, (27) * 


* 


% 


KKK KKK KKK KKK KK KKK KK KK KK KKK KEK KERR ERK K KEKE KKK KEKE KER KK KKK KEKE RK KERR KER ERE RR 
SUBROUTINE CDAT85( LOW, HIGH, VALUE, NAME, IMOD) 
COMMON/FILS/IIN, IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 
REAL LOW 
CHARACTER*35 NAME 
CHARACTER*1 ESC 
ESC = CHAR(27) 

10 FORMAT.” VENTER THE =A3o) 


20 FORMAT(’ THE MINIMUM VALUE FOR ’,A35,/' IS ',F6.2) 
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30 FORMAT(' THE MAXIMUM VALUE FOR ',A35,/' IS ‘',F6.2) 


FKKHK HK HHH HK AK KKK HK KH KEK HE KH KE KH KEK KEKE KKK KH KH BEAK KKK HHH HK KKH KK HH HK KKK KH KKH 


* * 
* CLEAR SCREEN, POSITION CURSOR, AND * 
* SET IMOD FLAG AND WRITE PROMPT TO SCREEN. * 
* * 


KKKKK KEK KH KH HH HK HK KEK KEK HK KEK HE KH KKH K HH KHK Hk KH Kk Wok Kk K ok Kok Kk Kok kK HK AK HK KKK KK 
100 WRITE (IOUT,*) ESC,’[2J' 

WRITE (IOUT,*) ESC,’[10;0H' 

IMOD = 1 

WRITE(IOUT,10) NAME 


KRKEKEKEKKEKKRKE KEK KEK KEKE KE KAKA E EEK KE KKK EE KE KARE KE EEA KK EEK KKK KEK KEKE LEK EEK SE 


* x 
* READ SCREEN VALUE. ON ERROR OR END OF DATA SET, GO BACK TO * 
* PROMPT. : 
* * 


KRKEKEAKKAEKEK KEKE KEK KEK KE KK HE KK KEK HE EK KEE KK EK KKK KEKE KEE KEK KEKE EEK H KKH HHH 


READ(IIN,*,ERR = 100, END = 100) VALUE 


KHKEKHKE KEKE BHA KKK KKH KH KKH HHH KKK AK E KEK HHH HH HEH HK HHH HHH HH HHH HHH HHH WAH HK HK AK 


* * 
* IF VALUES ARE OUT OF RANGE , PRINT ERROR MSG AND RETURN TO * 
* PROMPT. * 
* * 


KHHKKK KKK KH KKK KK KEKE KKK KEELE KKK KKK KKK HERB RAKHI BH HH HEHEHE HEHE HEHEHE HE HEHE 
IF (VALUE.LT.LOW) THEN 


WRITE (IOUT,20) NAME, LOW 


ye) 


CALL DELAY (6) 
GO TO 100 
ELSE 
IF(VALUE.GT.HIGH) THEN 
WRITE (IOUT,30) NAME,HIGH 
CALL DELAY (6) 
GO TO 100 
ELSE 
END IF 
END IF 
RETURN 


END 
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GAUSS 


KKK KKK KKK AKER AKAKKE KEKE KE KK KEKE KK KKK KK KEK KK KEKE KK EK KEKE KEKE ER ERK RER EE RERE 


* 


* 


* 


SUBROUTINE GAUSS (NRHS ) 


SOLUTION OF LINEAR ALGEGRAIC SYSTEM BY GAUSS ELIMINATION 
WITH PARTIAL PIVOTING. TAKEN FROM ‘AN INTRODUCTION TO 
THEORETICAL AND COMPUTATIONAL AERODYNAMICS’, JACK MORAN, 


JOHN WILEY & SONS, NEW YORK, 1984, PG 78. 


VARIABLES USED ARE: 
A() : COEFFICIENT MATRIX, AUGMENTED WITH RIGHT HAND 
SLDES IN COLUMN NEQNS+1, INPUT/OUTPUT, UNUSABLE 
AFTER RETURN 
NEQNS : NUMBER OF EQUATIONS IN MATRIX, INPUT 


NRHS : NUMBER OF RIGHT HAND SIDES, INPUT 


* 


* 


* 


* 


* 


* 


BHKKKK KEK AA KK KKK KKK KEKE AAKK KE KEKE KEKE KKK KEKE KKK KEK AK KEKE KK KEKE KK EEE KEKE HK 


SUBROUTINE GAUSS (NRHS) 
COMMON/COF/A(101,101) ,NEQNS 
NP = NEQNS+1 


NTOT = NEQNS+NRHS 


KEKKKK KEKE K KEKE KEK KEK KKK KEKE KEK KKK EKER EKER KEE RK KERR ERR EK KK ERK ERK ERK EEK EN 


* 


* 


* 


SEARCH FOR THE LARGEST ENTRY IN THE (1I-1)TH COLUMN, ON OR 


BELOW THE MAIN DIAGONAL. 
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x 


* 


* 


* x 
FH HAHA AA HHH HAHAH HH HAHAHAHA HAA HEE EE HEE HEHEHE DEDEDE EEA FEE FEE HEHE HH HH HHH 
DO@lS0 I = 2; NEONS 


IM = I-l 


: 


EM 


; 


ABS (A(CIM, IM) ) 
DO 110 J gee 1, NEOQNS 
IF (AMAX.GE.ABS(A(J,IM))) GO TO 110 
IMAX = J 
AMAX = ABS(A(J,IM)) 
LO CONTINUE 


KEKKEEKAKAK KEKE KK KEKE KK KEK KEKE KKK KKK KE EK KEK KEK EEK EK KEKE KEE KEKE KEKE KKK HK KEKE KKK EK 


* * 
* SWITCH THE (I-1)TH AND IMAXTH EQUATIONS * 
* x* 


HARK EKKKK KKK KK KKK KAKKKKEKE KK KK KK KK KKKKK Kk Kk Kk ok ok kk Kk Kk Kk KKK KK kw KK Kk Kk KK 
IF (IMAX.NE.IM) GO TO 140 
DO 130 J = IM,NTOT 
TEMP = A(IM,J) 
A(IM,J) = ACIMAX,J) 
A(IMAX,J) = TEMP 
130 CONTINUE 


KKKKK KKH K HHH AH HK HK AKA HK KKH HH HH HK HK HII HEE EK IEEE Oe IE KE EK EK ak ok ak Ok aE EE AK a HK OE HK HEA KN K 


* * 
* ELIMINATE THE (I-1)TH UNKNOWN FROM ITH THROUGH (NEQNS)TH * 
* EQUATIONS. * 
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og * 
PKK KKK KH RHHKK HH HHKHKK HKHKHKHH KH HH KK HK KEKE AK AK HAKKAR KAKA KEKE KER KEKE KEK 
140 DO 150 J = I,NEQNS 

R = A(J,IM)/ACIM, IM) 

DO 150 K = I,NTOT 

A(J,K) = A(J,K)-R*ACIM,K) 
150 CONTINUE 


KREKAAKKEKKAKKKEKEKERK KEKE KEREKEKKEERKKAEKRERE KERR RE R RE RR ERR REE RERREREEEN 


* * 
* BACK SUBSTITUTE * 
* * 


KREKEKKEKRA KEKE KKK KEKE KEK EE KKK RE KEKE EKER EKER KEK AKER KE KEKE KEK KK KE KEKE KEKE 
DO 220 K = NP,NTOT 
A(NEQNS ,K) = A(NEQNS ,K) /A(NEQNS , NEQNS) 
DO 210 L = 2,NEQNS 
I = NEQNS+1-L 
Je ilar 
DO 200 J = IP,NEQNS 
AC(I,K) = ACI,K) - ACI,J)*A(CJ,K) 
200 CONTINUE 
A(I,K) = A(I,K)/A(I,1) 
210 CONTINUE 
220 CONTINUE 
RETURN 


END 


79 


RNDUP 


KKH HK AHHH WK HAW WH HH MMH HHH HW HHH WH HW HK HH WKH HWM WH HW WH HH WKH HH HW Wok KH Hk ok Hk HW ak HK HK HA 


* * 
* SUBROUTINE TO ROUND UP THE VALUE OF AN AXIS LIMIT. * 
* FOR EXAMPLE IF THE VALUE PASSED IS .0018, THE VALUE RETURNED * 
* IS .0020. THIS IS NECESSARY TO ENSURE THAT AXIS LIMIT TSeArE * 
* WELL DEFINED. * 
* * 
* VARIABLES USED ARE: * 
* VALIN THE VALUE PASSED * 
* VALOUT THE VALUE RETURNED * 
* VAL A TEMPORARY VARIABLE * 
* IRTN AN INTEGER MULTIPLIER USED TO RETURN 1G THE CORRECT * 
* ORDER OF MAGNITUDE. * 
* ISIDE AN INTEGER ADDED TO ROUND UP. DEPENDS ON WHETHER THE * 
* VALUE PASSED IS POSITIVE OR NEGATIVE. * 
* cv. A TEMPORARY VARIABLE * 
* * 


HRKK HHH HHH HAHAHAHAHAHA HHH HHH HAHAHAHAHAHA AHHH HEE DEDEDE EE HE HHH HHH TBE BEBE DE EE HH EH 
SUBROUTINE RNDUP(VALIN, VALOUT) 
RIN = 1.0 


VAL 


VALIN 

IF (VAL.EQ.0.0) THEN 
VALOUT = VAL 
GOTO 200 


ELSE 
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100 


200 


ENDIF 


IF (VAL.LT.0.0) THEN 


ISIDE ane 


ELSE 


I 
me 


IESE BY D 
ENDIF 


IF (ABS(VAL) .GE.10.0) THEN 


VAL = VAL/10.0 
RTN = RTN*10.0 
ELSE 


Cabo (VAI). LT. 1.0) THEN 


VAL 


VAL*10.0 


RTN 


RTN/10.0 
ELSE 
IV = AINT(VAL)+ISIDE 
VALOUT = REAL(IV*RTN) 
GO TO 200 
ENDIF 
ENDIF 
GO TO 100 
RETURN 


END 
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SDIG 


HREKKKK KK KK EK KKK KKK KKK HE KKK KK KKK KKK HHH HHH EKER HEHEHE IEE EEE HEHE HEHE DE DEE DE EE IE HEE IE HEE EH 


* * 
* SUBROUTINE SDIG(MININ, MINOUT , MAXIN, MAXOUT, IDIG) * 
* * 
* SUBROUTINE TO GET THE RANGE OF VALUES IN THE AXIS LIMITS * 
* TO BE 5 UNITS APART. AND RETURN THE NUMBER OF SIGNIFICANT * 
* DIGITS FOR THE AXIS. * 
* * 
* VARIABLES USED ARE: * 
* MININ THE MINIMUM AXIS VALUE PASSED * 
* MINOUT THE MINIMUM AXIS VALUE RETURNED * 
* MAXIN THE MAXIMUM AXIS VALUE PASSED * 
* MAXOUT THE MAXIMUM AXIS VALUE RETURNED * 
* VAL A TEMPORARY VARIABLE * 
* VALNOT A TEMPORARY VARIABLE * 
* RTN A MULTIPLIER USED TO RETURN TO THE CORRECT * 
* ORDER OF MAGNITUDE. * 
* ISIDE AN INTEGER TO DETERMINE WHETHER THE VALUE IS POSITIVE * 
* OR NEGATIVE . 
* IV A TEMPORARY VARIABLE * 
x * 


HHH RHR HER HEH KHAN HAHAHAHAHA HHH HEHEHE HH HEH HHH EEF EDE EEE HH HH HH HEHE K WH 
SUBROUTINE SDIG(MININ,MINOUT , MAXIN ,MAXOUT, IDIG) 
REAL MININ,MINOUT , MAXIN,MAXOUT,RTN, VAL, VALNOT 


INTEGER IDiG 1CGASE TSlDE A eGo 
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RIN = 1.0 
IF(ABS(MININ) .LT.ABS(MAXIN)) THEN 


FF KHH AHHH HHH AHHH HHH AHHH HHH HHH HEH HHH HH HK HHH HHH HH BE sk ak HW Me se ak 9k se ak ak ak Hk vk WE ak Hk Hk kK HE I HK 


* * 
* ICASE = 1 IS FOR THE MIN AXIS VALUE BEING SMALLER IN MAGNITUDE * 
* THAN THE MAX AXIS VALUE. * 
* * 


HHAKHKAKKKKK KKK HEH K HHH KKK HR HHH HHH HHH HK HHH HH HK de Hk ke MH mk vk He ak vk Hk ak vk ok vk Hk ek He xk HK 


ICASE = l 


VAL = MAXIN 


VALNOT = MININ 


ELSE 


BKKKKKHKHKK KK HHH HK KHANH KKK KKHKHH HK HHH KK HK KKK KKK HK KKH AK KH HHH KKK KKK KKK KKK KK KKH 


* * 
* ICASE = 2 IS FOR THE MAX AXIS VALUE BEING SMALLER IN MAGNITUDE * 
* THAN THE MIN AXIS VALUE. * 
* * 


KHHAH AKA HA KAKA KKK KKH KK HH KH HEH HEHE HHH KK HHH HHH HE HHH AHHH HHH HK HK HH HAH HK  k ak ke vk kk Kk H 


ICASE = 2 


VAL = MININ 


VALNOT = MAXIN 


ENDIF 


FHHHHKHK KKH HHH HHH HHH HHH HHH HHH HHA HHH HHH HHH HHH HHH HHH HHH HHH HH HHH HHH HHH KKK KK 


* * 
* DETERMINE WHETHER THE SMALLER (IN MAGNITUDE) AXIS VALUE IS * 
* POSITIVE OR NEGATIVE. ss 
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¥ * 
KHREBKEKEKEBKE KAKA REE K KBE BK KKK KEKE KKK EKER EE eke kk ek hbk ok bk b&b ok ob ak ok ok 
ISIDE = SIGN(1,VALNOT) 
100 IF (ABS(VAL).GE.10.0) THEN 
VAL = VAL/10.0 
RTN = RTN*10.0 
ELSE 


LF (ABS (VAL een 1.0). THEN 


VAL = VAL*10.0 
RIN = RTN/10.0 
ELSE 
IF (VALNOT.EQ.0.0) GOTO 300 
IV = NINT(VALNOT/RTN) 
LF(IV.EQ.0) IV = IV + ISIDE 
GO TO 200 
ENDIF 
ENDIF 
GO TO 100 


200 CONTINUE 
ITEST = ABS(NINT(VAL) -IV) 
IF ((ITEST.EQ.5).OR. (ITEST.EQ.10).OR. (ITEST.EQ.15).OR. 


* (ITEST.EQ.20).OR. (NINT(VAL) .EQ.-IV)) GOTO 300 


IV = IV + ISIDE 
ITEST = ABS(NINT(VAL) -IV) 
IF (CITEST.EQ.5) .OR. CITEST. EQ,10) ORV CPTES Fehon eo 


* (ITEST.EQ.20) .OR. (NINT(VAL) .EQ.-IV)) GOTO 300 


1° 


VAL = VAL - ISIDE 
ITEST = ABS(NINT(VAL) -IV) 
IF ((ITEST.EQ.5).OR. (ITEST.EQ.10).OR. (ITEST.EQ.15).OR. 
* (ITEST.EQ.20) .OR. (NINT(VAL) .EQ.-IV)) GOTO 300 
IV = IV + ISIDE 
ITEST = ABS(NINT(VAL) - IV) 
IF ((ITEST.EQ.5).OR. (ITEST.EQ.10).OR. (ITEST.EQ.15).OR. 
* (ITEST.EQ.20).OR. (NINT(VAL) .EQ.-IV)) GOTO 300 
VAL = VAL - ISIDE 
300 IDIG = NINT(LOG10(1.0/RTN)) + 1 


ir CLCASE. EQ’. 1), EHEN 


MINOUT = IV*RTN 
MAXOUT = NINT(VAL)*RTN 
ELSE 
MAXOUT = IV*RTIN 
MINOUT = NINT(VAL)*RTN 
ENDIF 
RETURN 
END 
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AR 


CMac 


APPENDIX B. NOMENCLATURE 


aspect ratio 

span 

local enormd 

average chord 

root “eh ome 

tip chord 

drag coefficient 

inauced drag coerrticiene 
it tse oe Pic iene 
tdealllite coefirerenc 
about the 


moment coefficient 


leading edge 
aie fou 
root chord for circulation=imede 


moment coefficient 


moment coefficient about the y axis 


of 


moment coefficient about the aerodynamic cena 


pressure difference coefficient 
grid element inset 

drag 

section induced drag 

FOTrcem per unit Sspan 

Vortex adasag ~iaicwiom 


Lies 


8:6 


ait 


es y) 
(Xo 5Vo) 
a 

fF 


oe 


moment 


@dastancembetween field and control points 


planform area 

local effective velocity 
remote velocity 

induced velocity 
chordwise coordinate 


x coordinate of aerodynamic center 


chordwise position of center of pressure 


Spanwise coordinate 

height coordinate 

field point 

eonerol point 

angle of attack 

Cr ec wila econ 

imeremental circulation 

leading edge sweep angle 
transformed spanwise coordinate 
tangent of leading edge sweep angle 
model 

Sapper facto ,in VORTEX model 


density 


ita 


MeEecMemwine “caper! in circulation model 


Peanstormed chordwise coordinate 
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Clive lat on 


CO (tie O lm © le position where velocity is evaluated 


Lield  poume position of the vortex element 
wing shape surface slope of the wing, including twist 


from root to tip and camber 
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